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1,0 Introduction 



This report documents a program for the solution of algebraic and 
implicitly defined stiff differential equations. We were particularly 
interested in solving very large systems of differential equations 
arising from partial differential equations where the finite element 
method has been applied in the space variables. 

Our original goal was to use a compact storage scheme for the 
large matrices involved and to use iteration to solve the linear 
algebraic systems which occur. However, the resulting program is easily 
adapted to different applications through user modifications accomplished 
by replacement of one or two relatively simple subroutines. Thus the 
program is a powerful one which can be used in a variety of applications. 
Four examples, illustrating different matrix storage techniques and 
different linear equation solvers are given in the appendix. Other 
storage schemes and solution methods, e.g. symmetric Jacobian with 
Cholesky decomposition, are relatively simple to implement. 

In Section 2 a brief review of the integration scheme is given. 

In Section 3 a discussion of differences between this program and the 
one from which it was adapted is given. Section 4 is devoted to a 
discussion of information concerning the use of this package. 

2.0 Theoretical Background 

Consider the system of N = m + il differential and algebraic 
equations in Vi** * * 

(1) F(y,y,t) + P(t) V = 0 , 



1 



with all or some of the initial values y^^Ctg) , . . . ,y^(tQ) , . . . , 

v^(tg) specified. Enough of the above values must be given in order 

to determine the remaining values and initial values for any of the 

• • 

derivatives, appears in equation (1). In equation (1) , 

P(t) is an N X matrix, F is a vector with N components, and 
V is a vector. 

The program documented here is a modification of a program due to 
Brown and Gear [1]. The method of solution is a modification of Gear’s 
method for stiff differential equations [2]. The application to differ- 
ential algebraic systems was given by Gear [3]. We will briefly describe 
the method here for completeness, and refer the reader to the references 
for more details. 

Suppose that approximate solution values are known at a number of 
equally spaced points, ^n-l’^n 2****^^n-k* represented by 

y^’^ ^\...,y^^ respectively. Let ^(^^-1^ represented by 

y(n 1) ^ ^ backward differentiation formula gives 

, • (n) 1 f (n) (n-1) (n-k) v 

h y^ = (ttg y + y + . . . + otj^ y ) y 

where h = coefficients and 6 q are from Gear 

[2] , p. 217. Substitution of this into (1) gives 



(2) ^ , t^) + P(t^) = 0 



n 



as the equation which y^'^^ and must satisfy. In equation (2) , 



V 1 / (n-1) . . (n-k). 

+ ... + ai,y ) • 
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In general, equation (2) represents a system of nonlinear equations for 
and . The method used for solving this system of algebraic 

equations is a variant of Newton's method. The initial guess, , 

is obtained by polynomial extrapolation using a Hermite polynomial 
interpolating the known values y^^ y^'^ ..., y^^ . Thus 

the predicted values has the form 



(3) 



(n),0 , 

y = h 



6, y 



(n-l) 



+ a, y 



(n-l) 



a y 
n-k ■' 



(n-k) 



The application of Newton's Method to equation (2) then yields the 
corrector equation 



(4) J • 



(n) , i+1 _ (n) , i 

^(n) ,i+l _ ^(n) ,i / ^0 



where J is the Jacobian matrix. 



(5) 



ay 3y 




Gear shows that the initial guess for is not important, 

and thus is used. Up to three iterations are performed on the 

corrector equation. The matrix J is not evaluated at each iteration, 
nor even at each timestep. J is evaluated whenever (i) the timestep 
or order is changed, or (ii) the corrector iteration fails to converge 
in three iterations. 

If the corrector iteration fails to converge, the J matrix is 
evaluated, unless it had already been evaluated at the current time. 
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If it has been evaluated at the current time, the timestep is reduced 
by a factor of 4, In either case the step is then retired. 

If the corrector iteration converges, the local error is estimated, 
based on the fact that the local error is approximately proportional to 
the difference between the predicted and corrected values of y^^^ • 

For this purpose a relative error tolerance is used for large solutions 
and an absolute error tolerance for small solutions. The root-mean- 
square norm (Euclidean norm divided by the square root of the number of 
components) is used for the vector with components e^/ymax^ , where 
e^ is the estimated local error in ymax^ = max (|yf^^|, 1) . 

If the error is larger than that specified by the user, an acceptable 
timestep is estimated for the current order or order one lower, and 
the step repeated. Up to three such failures are permitted, after which 
an attempt is made to start over with a first order method. 

When using a method of order q , the program takes at least 
q + 1 steps before changing the timestep. Changes in timestep are 
preceeded by calculation of the predicted timestep at current order and 
order one higher and one lower. If the timestep can be increased by more 
than 10%, the order corresponding to the largest possible timestep is 
used. If the timestep cannot be increased by at least 10%, the current 
order and stepsize are retained for at least 10 more steps. 

After each step a test is made to determine whether time has 
advanced to or beyond the input end time. Control is returned to the 
calling program if it has. 
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At the initial call, no history of the solution is available, so 
the program must begin with a first order method, taking two such steps 
in accordance with the above description. The times tep must be suitably 
small, again in accordance with the above. At the point the program can 
begin to increase the order of the method and the times tep. Because 
the Jacobian matrix must be generated whenever the timestep is changed, 
it is not efficient to try too large a timestep initially. Because the 
program rapidly finds the best order and timestep, it is relatively cheap 
to underestimate the initial timestep compared to the cost of overestimating 
it. 



3,0 Differences compared with DFASUB 

The principal differences between the SDESOL/LDASUB package and 
DFASUB, and the reasons for incorporating them are as follows. 

(i) The main goal of this revision was to generate a program which 
could solve very large sparse systems of differential equations efficiently, 
both in terms of storage requirements and execution time. We are 
particularly interested in the solution of ordinary differential equa- 
tions arising from partial differential equations where the finite 
element method has been used to discretize the problem in space. 

Large sparse problems require at least a different system of stor- 
ing the Jacobian matrix and possibly the use of iteration to solve for 
the quasi-Newton iterates in equation (2.4). Two such subroutine 
packages, to be used with the basic subroutines, have been provided. 

Another package using standard elimination techniques is also provided 
and is convenient for smaller systems of equations. Use of any of these 
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options requires the user to supply a subroutine to evaluate his form 
of the equations (1), and for efficiency, a subroutine to explicitly 
evaluate the Jacobian • A subroutine to approximate a full Jacobian 
through numerical differencing has been provided. With the exception 
of a minor correction, this is the same as given in [1]. While use of 
this routine is convenient, it is inefficient and should be avoided 
for large systems. It is anticipated that the user can provide his own 
subroutine package, using his own storage scheme for the Jacobian, and 
with a suitable equation solver for the Newton iterates. There should 
be no need to disturb the basic package which carries out the time 
integration. 

(ii) For user convenience, without a major rewrite of DFASUB, a driver 
routine, SDESOL, to be referenced by the user and which then communicates 
with LDASUB was written. The chief function of SDESOL is to set up a 
number of references to work storage areas required by LDASUB. In 
addition, some testing of parameters is accomplished, and a subroutine 

to calculate initial values of derivatives is called. 

(iii) A subroutine to calculate initial values of derivatives, DERVAL, 

has been provided. The routine provided requires that the first m 
9F 

rows of — be nonsingular, which does not need to be true in the 

9y 

general case. For this reason, and others discussed in Section 4, the 
user may need to provide either his own version of DERVAL to evaluate 
the derivatives initially, or else he may supply initial values and a 
dummy version of DERVAL. 

(iv) Other changes made in generating LDASUB from DFASUB were to simplify 
the code for the particular type of problem we wish to solve, while 
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others were to enhance the usability of the code. Some errors were 
also corrected, notably two errors in coefficients for the fifth and 
sixth order methods. DFASUB had the capability of computing various 
elements of the Jacobian at different times if they had different 
dependencies, with the possibility of inverting that part of the matrix 
at that time, if it could be done. This could result in increased 
efficiency in certain problems, at some expense in convenience, but for 
our purposes it was not considered useful, and was removed. Therefore 
only one call is made to evaluate the Jacobian. The Jacobian was 
evaluated at the beginning of each times tep in DFASUB, but this has 
been eliminated in LDASUB, in accordance with the description in Section 
2. A subroutine, S2, was called from DFASUB to evaluate time dependent 
terms whenever time was changed. This is reasonable, since the function 
evaluation routine may be called several times at a given value of the 
independent variable. We have removed this, preferring to test for a 
new time in the routines where time dependent terms appear, then 
evaluating and storing them internally to that routine when necessary. 
This helps make the function evaluation more self contained, as well. 

In DFASUB extra parameters in the calling sequence allowed the 
user to communicate constants to the function evaluation and Jacobian 
subroutines via DFASUB. We believe this is inefficient and confusing, 
and removed this capability, preferring to communicate from the main 
program to these subroutines via Common, or possibly through multiple 
entry points. 

The norm used for error tests in DFASUB is the Euclidean norm. 

This has the undesirable property that for large systems the allowable 
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error criterion may be large. We therefore changed to the root^mean- 

square norm in LDASUB, which is simply the Euclidean norm divided by 

the square root of the number of components. One other change was 

made in the error tests. As noted in Section 2, the error vector has 

components e^/ymax^ , where e^ is the estimated local error in the 

i— variable y. , and ymax, = max (|y^^^|, 1) . In DFASUB the 
^ 0£k_<n ^ 

maximum was taken only up to the previous timestep, n-1 . This change 

was incorporated because some of the problems in which we were interested 

began with many components at zero, but which very rapidly became large, 

12 

around 10 or more. Without updating the value of ymax^ , the 
size of the timestep was artificially kept extremely small in order to 
satisfy an unreasonable error tolerance. For this reason, the maximum 
value of the component was updated before the norm of the relative errors 
was computed. For problems where values range near to one, the modifica- 
tion will result in no appreciable change in performance. 

The printout of counters, timestep, time, and values of the dependent 
variables was made an option through a parameter in the calling sequence. 

An additional value printed is the order of the method being used. 

DFASUB incorporated the capability of terminating if a certain 
number of floating point overflows had occurred. This capability was 
removed from LDASUB. 

The final modification to the program was the incorporation of a 
restart capability without having to begin again with a first order 
method. This was accomplished by adding two entry points to LDASUB. 

One, LDASAV, saves values internal to LDASUB, returning them to the 

main program,- where they can be saved for the time at which the calculation 
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is to be resumed. At that time, another entry point, LDARST, restores 
those values internal to LDASUB, while other necessary values are 
restored through the calling sequence. 

4.0 Subroutine Descriptions 

The description of subroutines is divided into two subsections. 

The first deals with the basic integration routine and other subroutines 
which make up the core package, and which should not be modified by the 
casual user. The second deals with a set of supporting subroutines, at 
least one of which must be supplied by the user since it defines his 
system of equations. The others may be usable in the form given in 
one of the examples, or can be defined by the user to accomplish his 
desired implementation. 

4. 1 Basic Subroutine Package 

4.1.1 Subroutine SDESQL 

This routine is the only one which needs to be referenced by most 
users. It serves as a driver for the integration routine, LDASUB. 

SDESOL has a simpler calling sequence than LDASUB and relieves the 
user of having to set up a number of auxiliary storage arrays. In 
addition, the routine calls DERVAL to calculate the values of the 
derivatives on the initial call. 

The calling sequence is 

CALL SDESOL (Y,YL,T, TEND, NY, NL,M,JSKF,MAXDER,IPRT,H,HMIN,HMAX,RMSEPS,W) 
where the parameters are defined as follows. 
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Y - Input and output. An array dimensioned (7, NY). On the initial 
call this array contains the initial values of the dependent 



stored in Y(j+l,i) , Here h is the current stepsize. 

These values must not be changed between returns to the 
calling program and subsequent entries to SDESOL. It is 
possible to interpolate for values of the dependent variables 
at times other than those calculated by using the formula 



currently in use, and is obtained as q = |JSKF/10| . 

YL - input and output. Array of linear variables, v^, i=l> •••, ^ 
The user supplies initial values for these variables, and 
during execution it contains current values of the linear 
variables. 

T - input and output. The user supplies the initial time, which 
is updated to current time during execution. 

TEND -input* Time at which the integration is to end. This is the 
only parameter normally changed by the user between succes- 
sive entries to SDESOL. 

NY - input. The number of differential and nonlinear variables, m 

NL - input. The number of linear variables, H . This may be zero. 

M - input. The number of variables to be included in the local 

error test. The error test will be performed for the variables 
y^, i=l, M . The M used is no greater than NY . 





j=0 



, where q is the order formula 
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JSKF -input and output. An indicator: on input, 

JSKF = 0 indicates that this is the initial call to SDESOL. 
Initial values of the derivatives are calculated and auxiliary 
storage references are set up. This also indicates to 
subroutine LDASUB to initialize parameters and begin with 
a first order integration method. 

JSKF > 0 indicates a continuation of a previous call to SDESOL 
JSKF = - 1 indicates a restart call. This is discussed 
further in Section 4.1.2. 

JSKF < - 1 may result from the user neglecting to test for 
error returns from SDESOL. Because of this possibility, 
the run is terminated with an appropriate comment when 
JSKF < “ 1 is input. 

On output, JSKF normally is a two digit number, ± qp . q 
is the order of the formula currently being used for the 
integration, p is an indicator determining the type of 
return. JSKF > 0 , p = 1 is the normal return. Note that 
SDESOL may be re-entered to continue the solution without 
changing JSKF. JSKF <0 is an error return, with the value 
of p indicating the error, as follows, 
p = 1 error test failure for H ^ HMIN 
p = 3 corrector failed to converge for H > HMIN 

p = 4 corrector failed to converge for a first order method 

p = 5 error return from subroutine NUITSL 

p = 6 error return from subroutine DERVAL 
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MAXDER- input . Maximum order method which should be used. The 
highest order possible is six. 

IPRT- input. Print control indicator. 

0 , no print from LDASUB 

> 0 , at each step, print number of steps, number of Jacobian 
evaluations , current order being used, stepsize for next step, 
current time, and current values of the dependent variables. 

H - input and output. On initial call it is an estimate of the 

timestep. During execution it is updated to the current value, 
and on return contains the stepsize to be tried for the next 
step. The input value need not be accurate. It is better to 
underestimate than to overestimate the initial value. The 
stepsize and order are varied to meet the local error tolerance 
specified. The user does not normally change the stepsize 
between entries to SDESOL. 

HMIN- input. The minimum stepsize to be allowed. 

HMAX- input. The maximum stepsize to be allowed. 

RMSEPS-input . The error test constant. The values of the relative 
local errors must have root -mean-square norm less than RMSEPS. 

W - an array of auxiliary storage required by LDASUB. This array 
must contain a total number of locations equal to the sum of 
(i) 13*NY + 5*NL for arrays used in LDASUB, (ii) storage 

for. the Jacobian matrix, and (iii) any locations used in 
processing the Jacobian, e.g., scratch storage used by an 
equation solver. 
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4. 1,2 Subroutine LDASUB 



This subroutine is the basic integration routine and performs the 
process in essentially the same manner as subroutine DFASUB. A brief 
description is given in Section 2 and differences between this routine 
and DFASUB are outlined in Section 3. Parameters in this routine in 
which the user may be interested are stored in the W array, an argument 
of subroutine SDESOL, 

YMAX - array of maximum magnitudes of the independent variables, 
y^, up to the current time (or one, if less than one). 

This is stored beginning at location 7*NY + NL + 1 of 
the W array. 

ER “ This is the array of differences between the predicted 
and corrected values of the variables, y^ , and is 
proportional to the estimated error. This array is 
stored beginning in location 8*NY + NL + 1 of the W 
array. 

This subroutine incorporates a restart capability. In order to 
restart from a previous point without beginning again with a first order 
method, it is necessary to have saved a number of variables internal to 
LDASUB, and then restore them before calling SDESOL again. To save the 
internal parameters, the user calls subroutine LDASAV(SAV). Here SAV 
is an array of length 29 in which the values to be saved will be stored. 
In addition to SAV, the user must also save a number of arrays in the 
calling sequence of LDASUB, and this is most easily accomplished by 
saving the W array in the calling sequence for SDESOL. Once these 
arrays have been saved, along with the other simple parameters in the 
calling sequence ( Y and YL need not be saved) , the user is free to 
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use the package to solve a different problem, or to terminate the computer 
run, to be restarted later. 

At the time the problem is to be restarted, the user calls sub- 
routine LDARST(SAV), where SAV is the array of values obtained previously 
by calling LDASAV. This restores internal values in LDASUB. The user 
then calls SDESOL with the same simple parameters and the W array as 
before, except that JSKF = - 1 and a new end time, TEND, is provided. 
Restoration of values (including Y and YL ) in LDASUB is completed 
and solution of the problem resumes. 

If the user desires to change the error tolerance, number of 
variables in the error test, or maximum order to be used, the user 
must make a new initial call to SDESOL, that is, set JSKF = 0 . 

4.1.3 Subroutine CQPYZ 

This subroutine simply transfers the contents of one array into 
another array. 

4.2 Supporting Subroutines 

This group of subroutinesmust , at least in part, be supplied by 
the user. The user must supply at least one subroutine, DIFFUN. For 
better efficiency, the user should supply a subroutine, JACMAT, to 
explicitly evaluate the Jacobian, although a version which approximates 
the Jacobian by numerical differencing is given in the appendix. To 
take advantage of sparsity or other features of his problem, the user 
will need to supply the subroutine NUITSL to solve the systems of 
equations (2.4). For certain problems the user may have to supply 
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subroutine DERVAL to calculate the initial values of the derivatives. 
We discuss the requirements of these subroutines in turn. 



4.2.1 Subroutine DIFFUN 

This subroutine simply evaluates the equations (2.1) at a given 
time and yalues of y , y , and V . Other parameters in the function 
definition must be transmitted from the calling program via COMMON or 
some other device, determined by the user. 

The calling sequence is 

CALL DIFFUN (Y, YL, T, HINV, DY) , where the parameters are defined as 
follows . 



Y 



YL 



T 

HINV 

DY 



input. Same as in SDESOL. This array contains the current 
values of the variables y^ and their (scaled) derivatives, 
input. Same as in SDESOL. This array contains the current 
values of the linear variables, 
input. Current time. 

input. 1/h , where h is the current stepsize. 
output. Array of function values. 



4.2.2 Subroutine JACMAT 

This subroutine evaluates the Jacobian matrix J , equation (2.5) 
at the given time and current values of the dependent variables, order, 
and stepsize. A version of JACMAT which approximates J by numerical 
differencing is given in the appendix. For maximum efficiency, the 
user should supply the explicit representation of the Jacobian. Because 
the Jacobian is used to solve for the quasi-Newton iterates, it is not 
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necessary for the Jacobian to be exact. Thus the user should consider 



the possibility of approximations which reduce the total number of 
computations in this step, with due regard for the fact that a smaller 
timestep may be required to obtain convergence of the corrector within 
three iterations. 

The calling sequence for this subroutine is 
CALL JACMAT (Y, YL, T, HINV, A2, N, NY, EPS, DY, FI, PW) , where the 
parameters are defined as follows. 



Y 


input. Same as in SDESOL, Y contains the current values 


YL 


of the variables y^ and their (scaled) derivatives, 
input. Same as in SDESOL. This array contains the current 
values of the linear variables. 


T 


input. Current time. 


HINV - 


input. 1/h , where h is the current stepsize. 


A2 


input. The constant from LDASUB. 


N 


input. Total number of variables. 


NY 


input. Number of differential and nonlinear variables. 


EPS 


input. Error constant from LDASUB, *RMSEPS. 


DY 


input. Array of current function values. 


FI 


scratch array of N locations available for use by this 
routine. 


PW 


output. The Jacobian matrix J , or an approximation, 
calculated in JACMAT and returned to calling program. 
This matrix is used in subroutine NUITSL and the storage 
mode must agree between the two subroutines. 
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4.2.3 Subroutine NUITSL 



This subroutine solves the equations (2.4) for the quasi-Newton 
iterates. This subroutine will normally be supplied by the user, 
although versions which solve the system by elimination methods and 
iterative methods, respectively, are given in the examples in the 
appendix. This subroutine will often be modified or replaced by the 
user to take advantage of sparsity or other features of his problem in 
connection with JACMAT, of course. 

The calling sequence for this subroutine is 
CALL NUITSL (PW, DY, FI, N, NY. EPS, YMAX, NEWPW, KRET) , where the 
parameters are defined as follows. 

PW - input. The Jacobian matrix J computed in JACMAT. 

DY “ input. Right hand side of the linear system to be 

solved . 

FI - output. The solution is returned in the array FI . 

N - input. Total number of variables. 

NY - input. Number of differential and nonlinear variables. 

EPS - input. Error constant from LDASUB, • RMSEPS . 

YMAX - input. Array of maximum magnitudes of y^ up to the 

current time (or one if maximum magnitude is less than one). 

NEWPW - input. Indicates whether a new J matrix has been 
computed since the last entry to NUITSL. 

- 1, indicates this is a new J matrix. If any 
preprocessing, such as LU decomposition is to be done, 
the preprocessing should be done and NEWPW set to zero. 

= 0, indicates the J matrix is the same as on the previous 
entry to NUITSL. 
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KRET - output. Return indicator. 

= 0 , normal return 

= 1 , error return, solution of equations not obtained. 

Note that the parameters EPS and YMAX are useful if an iterative 

method is used for solutuon of the equations. Because the solution 

represents corrections to the predicted value, and corrections to that, 

the solution is small compared to the dependent variable values. Hence, 

compared to the YMAX array, the error tolerance can be fairly large. 

The following convergence criteria have been used, with great success. 

th 

Let 6u^ denote the i — component of the difference between successive 

th 

iterates, with u^ being the i — component of the current iterate. 

Then the iteration is considered to have converged whenever 




Condition (i) requires convergence to 2 digits more accuracy than the 
user has asked for in the solution of the system (2.1), relative to 
YMAX. Condition (ii) requires the same relative accuracy in u^ as is 
asked for by the user in the solution of the system (2.1), unless the 
solution is smaller than c , in which case the change is compared to 
c rather than |u^| . This avoids difficulty if u^ is close to zero. 
The e above is EPS = v^*RMSEPS, where M and RMSEPS are inputs to 
SDESOL. Two versions of NUITSL incorporating iteration and this conver- 
gence test are given in the appendix. 
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4.2.4 Subroutine DERVAL 



This subroutine solves for, or otherwise supplies the initial values 

of the derivatives, and possibly other variables. In some instances 

it may need to be supplied by the user. The standard version of 

DERVAL given in the appendix uses Newton’s method to solve the first 

m (=NY) of the equations (2.1) for * assuming values for y(t^) 

8F 

and V(t^) have been supplied. To accomplish this, the matrix — 

-20 ^ 

is needed, and this is obtained by calling JACMAT with h = 16 , 

A2 = -1 , and N = NY , implying NL = 0 for this call. Special care 



must be taken if in fact NL is not zero to assure that the matrix is 

20 9F 

computed and stored properly. The matrix returned is then 16 — . 

9y 



A call to DIFFUN yields the function values 

F(y(to) , yC^q) 9 ^q) where yC^q) is the current iterate. 

20 

Multiplication of the function values by 16 and a call to NUITSL 

(again with N = NY) gives the Newton iterate. Of course, the same 

sort of special care as necessary in JACMAT is necessary in NUITSL. 

3F 

Obviously the above scheme cannot work if — is singular, such 

9y 

as it would be if one of the equations is algebraic. In this instance 
the user must either devise his own version of DERVAL, or supply the 
values along with a dummy version of DERVAL. In an extreme case the 
user may simply set initial derivatives to zero. This will provide a 
poor predicted value on the first step, and will force an artificially 
small times tep for the first two steps. However, the overall penalty 
is generally small, as appropriate (corrected) values are computed at 
the first step, and after two steps the program quickly increases the 
times tep. 
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The calling sequence for this subroutine is 
CALL DERVAL (Y, YL, T, N, NY, DY, KERET) , where the parameters are 
defined as follows. 

Y - input and output. Same as in SDESOL. On entry Y(l,i) 

contains the initial values of the variables . On 

return, the values of the derivatives are stored in Y(2 i) 

YL - input. Same as in SDESOL. This array contains the initial 
values of the linear variables. 

T - input, initial time 

N - input. Total number of variables. 

NY - input. Number of differential and nonlinear variables. 

W - The scratch array from SDESOL, can be used in any way 

needed by this subroutine. 

KERET - output. Return indicator 

= 0 normal return 

= 1 error return, initial values were not obtained. 
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Appendix 1: Program Listings 



The following are listings of the basic subroutine package and 
supporting subroutines which are of general use. For simple problems 
the user only needs to supply a calling program and a subroutine, DIFFUN, 
to evaluate the equations. Use of the NUITSL routine in computer facilities 
which do not subscribe to the IMSL package will necessitate modifications 
to replace LUDATF with another LU decomposition routine, and LUELMF 
with another forward and backward substitution routine. 
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ooooooooooooooooooooooooooooooooooooooooooooooonooooooooooooooooooooooooo 



SLBROUTINE SOESOL < Y , YL iT ,T END, NY ,NL , M, JS K F , M4X0E R , I PRT , H , HMI N 
IHRAX,RMSEPS,M 



CALL 



SLEPQUTINE SCESOL !S A DRIVER ROUTINE FOR SUBROUTINE LDASU6, 
ITS PURPOSE IS TO SET UP THE NECESSARY REFERENCES TO A LARGE 
BLOCK OF AUXILLARY STORAGE, AND OBTAIN INITIAL VALUES OF 
DERIVATIVES. 

THE CALLING SEQUENCE FOR SDESOL IS 



WHERE TFE PARAMETERS ARE DEFINED AS FOLLOWS. 

Y - ARRAY dimensioned (7,NY>. THIS ARRAY CONTAINS THE 

DEPENDENT VARIABLES AND THPIR SCALED DERIVATIVES. 
Y<J+l,n CONTAINS THE J-TH DERIVATIVE OF THE I-TH VA 
lABLE TIMES H*+J/ J-FACTGRI AL , WHERE H IS THE CURRENT 
STEPSIZE. ON FIRST ENTRY THE CALLER SUPPLIES THE 
INITIAL VALUES OF EACH VARIABLE IN YU, I). ON SUB- 
SEQUENT ENTRIES IT IS ASSUMED THE ARRAY HAS NOT 
BEEN CHANGED. TO INTERPOLATE TC NON-MESH POINTS, 
THESE VALUES CAN BE USED AS FCLLOWS. IF H IS THE 
CURRENT STEPSIZE AND VALUES AT TIME T+E ARE 
NEEDED, LET S = E /H AND THEN 



I-TH VARIABLE AT T+E IS 



JS 

SUM Y( J + l ,I )*S+*J 
J = 0 



YL 

T 

TEND 

NY 

NL 

M 

JSKF 



The value of JS IS OBTAINED IN THE CALLING PPCGRAN 
BY JS = IABS( JSKF/101 

ARRAY CF NL VARIABLES WHICH APPEAR LINEARLY. 

CURRENT VALUE OF THE INCEPENCFNT VARIABLE (TImcj 
END TIME 

NUMBER OF DIFFERENTIAL EQUATICNS AND NONLINEAR 
VARIABLES. 

NUMBER OF LINEAR VARIABLES 

NUMBER OF VARIABLES INCLUDED IN THE ERROR TEST 
AN INDICATOR USED 8CTH ON INPUT AND OuTPL'^ 

CN INPUT, JSKF = -I INDICATES A RESTART CALL TO 
SDESCL. JSKF * 0 INDICATES AN' INITIAL CALL TO 
SDESCL. JSKF > 0 INDICATES A CONTINUATICN OF THE 
PREVIOUS CALL TO SDESOL. JSKF < -I MAY HAVE RESULTS 
FROM THE USER NEGLECTING TO TEST FOR ERRCR RETURNS 
FROM SDESCL. BECAUSE OF THIS POSSIBILITY, JSKF < -1 
RESULTS IN TERMINATION CF THE RUN WITH THE 
APPROPRIATE COMMENT. 

ON OUTPUT, JSKF CONSISTS OF TwO DIGITS AND SIGN, 

♦OR - QP. Q IS THE ORDER OF THE FORMULA CURRENTLY 
BEING USED. ? INDICATES THE TYPE OF PETLRN, AS 



FOLLOWS. 
JSKF > 0, F 
JSKF < 0 IS 
MEANINGS, 



= I IS THE NORMAL RETURN 

AN ERROR RETURN, HITH THE FOLLOWING 



P = I ERRCR TEST FAILURE FOR H > HMN 

P = 3 CORRECTOR FAILED TC CONVERGE FOR H > HM! 

P = 4 CORRECTOR FAILED TC CONVERGE FOR FIRST 

ORDER METHOD 

P = 5 ERROR RETURN FROM SUBROUTINE >llirrsL 

P = 6 ERROR RETURN FROM SUBROUTINE OERVAL 

MAXDER - maximum ORDER DERIVATIVE THAT SHOULD BE USED IN 
METHOD. IT MUST BE NO GREATER THAN SIX. 

IPRT - INTERNAL PR’NT CONTROL INDICATOR FOR LDASU3. 

IPRT = 0 NO PRINT 

IPRT > 0 PRINT COUNTERS, STEPSIZE, CURRENT TI(' 
AND VALUES OF DEPENDENT VARIABLES AT 
EACH STEP. 

H - CURRENT STEPSIZE. AN INITIAL VALUE MlST BE SUPPLIEl 

BUT NEED NOT BE THE ONE WHICH MUST BE USED, SINCE T^ 
SUBROUTINE WILL CHOSE A SMALLER ONE IF NECCESSAPY TC 
KEEP THE ERROR PER STEP SMALLER THAN THE SPECIFIED 
VALUE. IT IS BETTER TO UNDERESTIMATE THE INITIAL 
STEPSIZE THAN TO OVERESTIMATE IT. THE STEPSIZE IS 
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SDE 


230 


SCE 


240 
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SDE 


29 C 


SCE 


200 


SDE 


210 


SCE 


220 


SCE 


230 


SOE 


34C 


SCE 


250 


SDE 


260 


SCE 


270 


SDE 


280 


SDE 


29C 


SCE 


400 


SDE 


410 


SDE 


420 
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SDE 


630 
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22 



ooooooooo on non ooooooooooooo on oooooooooooooooo 



NORI^ALLY NOT CHANGED PY THc LStR. SDE 760 

HMIN - MINIf^UM STEPSIZE ALLOWED 770 

H.MAX - MAXIMUM STEPSIZE ALLOWED -DE 780 

RMSEPS - THE ERROR TEST CONSTANT. THE PCjCT-m«= AN- SQUAR E OF SDE 7S0 

THE SINGLE STEP ERROR ESTIMATES, ER ( I ) , OTVID«=D 3Y 800 

YMAX(I) = ('MAXIMUM TO CURRENT TIME OF Yfll ) MUST 3E SDE 81C 

LESS THAN EPS. THE STcPSIZE AND/OR THE CRDEP SDE 82C 

ARE VARIED TO ACHIEVE THIS. SDE 830 

W - SCRATCH STORAGE ARRAY. MUST EE AT LEAST 13*N!Y ♦ 5»NLSDE 8^0 

LOCATIONS, PLUS THOSE REQUIRED FOR STORAGE OF THE SDE 850 

MATRIX PW (SEE DESCRIPTION CF SUBROUTINE JACMAT). SDE 860 

THE STORAGE JF PW WILL NORMALLY REQUIRE NO MORE THAN SCE 87C 

* 2^N LOCATIONS, AND IF COMPACT STnpAGF TECH- SCE 830 

NIQUES AP«= USED, CAN BE MUCH FEWER. SDE 890 

SCE 900 
-SDE 910 
:0t 920 

SDE 920 
SOE 9A0 
SCE 950 
SCE 960 
SDE 97C 
SDE 960 
SDE 990 

soe 1000 



DIMENSICN Y(7,l), YL(1), W(l) 

IF (JSKF.GT.O) GO TO iz6 
IF (JSKF.LT.-1) GC TD 140 
N = NY + N'L 

IF (JSKF.LT.O) GO TO 110 

IF THIS IS THE FIRST ENTRY, OBTAIN VALUES CF THE DERIVATIVES. 
CALL DFRVAL ( Y, YL , T ,N, NY . W , KRET R ) 

IF (KRETF.NE.O) GC TC 130 



THE ARRAY 
THE ARRAY 
THE ARRAY 
THE ARRAY 
THE ARRAY 
THE ARRAY 
THE ARRAY 
THE MATRIX 

110 NSVL = 7>>NY + 1 
NYMAX = NSVL^NL 
NER = NYMAX+NY 
NESV = .\EP + NY 
NFl = NESV+NY 
NCY = NFl+N 
NPW = NDY+N 
120 JS = JSKF 
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PW 


STAR TS 


AT 
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w 
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S C£ 
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CALL LDASUe (Y, YL , T , T FND * N , NY , ^ • JS , KF • M AX DER , I PRT , H , ! N ♦ HM A X , 

IRMSEPS, W, W<NSVL ) , w(NYMAX) , W ( Nr R ), W ( NESV ), W < N FI ), W ( NCY ), W ( NPW ) ) 



CODE JSKF ON RETURN FROM LOASUB 

JSKF ^ ISIGN(JS=<‘10+I ABS(KF )tKF) 
RETURN 

130 JSKF = -6 
RETURN 

140 PPINT 1, JSKF 
STCP 



1 FORMAT (*OIT IS AN ERROR T-0 ENTER SDES; 
1 » RUN HAS BEEN TERMINATED.') 

END 



•L WITH JSKF = ' , IlD// 
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1320 

1230 

1240 

1350 

126C 
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SLBROUTINE LCASUb ( Y , YL , T , TEND ,M , NY , M t JST ART , KFLA G , M A XDR , I PRT,Ht 
IHM IN, HM AX, RMSEPS, SAVE, YL3V ,YMAX, ER, cSV, FI, DY,PW) 

SIBRQUTIK£ LCASUb IS A MODIFICATION OF SUBRGUTINF DFASUB 

WHICH IS DUE TO R. L. BROWN AND C. W. GEAR. D'^ASUH IS DOCUMENTED 

IN THE REPORT 

C OCUMENTATIfN FOR 0FA5UB-- 
BY R. L. EROWN AND C. W. GEAR 
REPORT UILCDCS-R-73-575, JUlY 1973 
UNIVERSITY OF ILLINOIS AT UR EANA -C HAMP A I 3N 
UR3ANA, ILLINOIS 61801 
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LDA 

LCA 

LCA 



10 

20 

20 

40 

50 

60 

70 

80 

90 

100 

110 



23 



OOOOOOOOOOOOOOOOOOOOO ->00000000000000000000000000000000000000000000000000000 



THIS REPCRT IS AVAILABLE THE .NATIONAL TECHNICAL INFORMATION 

SERVICE CF THE U. S. DEPARTMENT OF COMMERCE UNDER ACCESSION NUM6E 
CC0-IA69-225. 

THE MODIFICATION HERE IS DOCUMENTED IN THE REPORT 

A PRCGRAM FOR THE NUMERICAL SOLUTION OF LARGE SPARSE SYSTEMS : 

ALGEBRAIC AND IMPLICITLY DEFINED STIFF DIFFERENTIAL EQUATIONS 

BY RICHARC FRANKE 

REPORT NPS53FE76051t MAY 1976 

NAVAL POSTGRADUATE SCHOOL 

MONTEREY, CALIFORNIA 939R0 



THE CALLING SEQUENCE FOR LDASUE IS 

CALL LDASUB (Y,YL, T,TEND,l^)>i'y,M.J.STAPT,KFLAG,MAXOR, IPRT,H,HMINt 
hMAX,RMSFPS ,S AVE , YLS V t YMA)C, E R , E$ V , F 1 , DY , P W ) 

k^HERE THE PARAMETERS ARE DEFI.NED AS FOLLOvJS* 

Y - ARRAY DIMENSIONED (7, NY). THIS ARRAY CONTAINS THE 

DEPENDENT VARIABLES AMD TJIEIR SCALED DERIVATIVES. 
Y(J+1,II CONTAINS THE j- fH DERIVATIVE OF THE I-TH VA 
IABlE TIMES H^*J/J-FAC^ORIALt *^HERE K IS THE CURRENT 
STEPSIZE. ?N FIRST ENTRY THE CALLER SUPPLIES THE 
INITIAL VALUES OF EACH VARIABLE IN Y(1,I) AwD AN 
ESTIMATE OF THE InIITIAL VALUES CF THE DERIVATIVES 
IM Y<2,I). ON SUBSEQUENT ENTRIES IT IS ASSUMED THAT 
THE ARRAY HAS NOT BEEN CHANGED. TO INTERPOLATE T^ 
NCN-MESH POINTS, THESE VALUES CA.N BE USED AS FOLLOWS 
IF H IS THE CURRENT STEPSIZE AND VALUES AT TIME T+E 
NEEDED, LET S = E /H AND THEN 



I-TH VARIABLE AT T+E IS 



NQ 

SUM Y( J+l ,I )*S*«J 
J = 0 



THE VALUE OF NQ IS 
BY NQ = JSTART. 



OBTAINED IM THE CALLING PROGRAM 



YL 



T 

TEND 

N 

NY 

M 



<0 



ARRAY OF NL = N - NY VARIABLES WHICH APPEAR LINEARLY 
THE USER SUPPLIES INITIAL VALUES FOR THESE VARIABLES 
CURRENT VALUE OF THE INDEPENDENT VAQZABLF (TIME) 

END TIME 

TOTAL NUMBER OF VARIABLES 

NUMBER CF DIFFERENTIAL EQUATIONS AND NONLINEAR 
VARI ABLES. 

NUMBER OF VARIABLES INCLUDED IN THE ERROR TEST. 

THIS NUMBER CAN BE NO GREATER THAN NY. IF IT IS 
GREATER THAN NV, NY VARIABLES ARE USED IN ERROR 

TFST. 

JSTART - INPUT AND OUTPUT INDICATOR. 

ON INPUT JSTART HAS THE FOLLOWING MEANINGS. 

THIS indicates a RE-START FROM A PREVIOUS 
POINT FOLLOWING TERMINATION OF THF PUN CR 
SOLUTION OF ANOTHER PRCBLFM DURING 1 HE SAME 
RUN. PARAMETERS IN THE CALLING SEQUENCE 
MUST HAVP been preserved FROM THE PREVIOUS 
USE, PARTICULARLY THE ARRAYS 
SAVE, YLSV, Esv, AND PW. 

THESE ARRAYS MUST BE SAVED AFTER 
TO SUBROUTINE LDASAV, WHICH ALSC 
NECESSARY PARAMETERS INTERNAL TC 
INDICATES AN INITIAL CALL TO LDASU6. THE 
ROUTINE INITIALI^FS ITSELF, SCALES thF 
DERIVATIVES IN Y(2,I) AND THEN PERFORMS THE 
INTEGRATION UNTIL T > TEND. 

INDICATES THE SOLUTION IS TO BE CONTIiNUED. 
AFTER THE INITIAL ENTRY IT lo NEITHER 
DESIRABLE NOR NCCESSARY TO RE-ENTER WITH 
JSTART = 0, SINCfc THIS RE- I NI T I AL ! ZES 
THE CODE, BcGMNFNG WITH A FIRST CRDER 
METHOD AGAIN. 

ON OUTPUT, JSTART IS SET TO THE VALUE OF NC, THE 
ORDER OF THE FORMULA CURRENTLY BEING USED. 
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or) 00 ooooooor> ooooooo ooooooooooooooooooooooooooooooooor»oor>ooor>oooono 



KFLAG 



MAX?*5 

IP^T 



THE CQNiPLETIGN 
MEANINGS 
+ 1 
-1 
-3 
-4 



:00c TNDICATU^, WITH THE POLLUTING 



THE INTEGRATION WAS 5LCCESSFUL 
cRROR TEST FAILURE FCR H > HMN 
CORRECTOR FAILED TO CCNVEPGE FO- H > 

CORRECTOR FAIlFC to CONVERGE FQ'^ FIRST 
ORDER METHOD 

-5 ERROR RETURN F‘=‘CM S'JBkCUTIME NUITSL 
MAXIMUM ORDER DERIVATIVE THA^ SHOUuD BE LSED IN THE 
method. it must be mo greater than Six. if IT TS 
GRFATF.R THAN SIX, THE MAXIMUM CRDEP U'ED NIlL BE SIX.LDA 



IC A 
I. CA 
LCA 
LDA 
LCA 
lCA 
LCA 
LCA 
LCA 
LCA 



HMIN 
HMAX 
RMS EPS 



SAVE 

YLSV 

YMAX 



INTERNAL PRINT CONTROL INPTCATCr 
^0 NO PRINT 

> 0 PRINT COUNTERS, '7EFSI2E, CURRENT -^-TME 
AND VALUES OF OE*'END£NT VARIABLES AT 

EACH STEP. 

CURRENT STEP5IZE. AN iNITiAL VALUE MUST BE SUPPLIED LDA 
BUT NEED MO-^ BE THE 1\S wH CH ^«ILL 6R USED. SIi^CE ThELCA 
SUBROUTINE WILL CHCCSE A SMALLER ONE :F NPCCESSARY TCL DA 
KEEP the error per STEP SMALLER THAN THE SPECIFIEC LCA 
value. it TS eETTEP TO UNDERF ST I MAT E THE INITIAL 
S-^EPSIZE than TC OVFkE:iTiMATE IT. THE STEPSIZE IS 
NORMALLY NOT CHANGED BY THE USER. 

MINIMUM STEPSIZE ALLOWED 
MAXIMUM STEPSIZE ALLOWED 

THE ERROR TEST CONSTANT. THE POOT-ME AN- SQUA ^ E OF 
THE SINGLE STEF ERROR ESTIMATES. £R ( I ) , DIVIDED EY 
YMAX(I) == <MAXIMUM TO CURR'=NT liME OF Y(T)) MUST PE 
LESS THAN RMSEPS. THE STEOSIZF AND/OR CHDEF APE 
VARI EO TO ACHIEVE THIS. 

AN ARRAY 0^ LENGTH AT LEAST 7*NY 
AN 4RRAY CF LENGTH AT LEAS"^ NL 

A VECTOR OF LFNGTH NY WHICH CONTAINS THE MAXTMyf/ 



LCA 

LCA 

LDA 

LCA 

LDA 



LCA 

LCA 

LDA 

LCA 

LDA 

LCA 

LDA 

LCA 

LCA 

LDA 

LCA 

LCA 

LDA 



OF EACH Y SE*=N SO FAR. CN THE FIRST CALL, THESF WiLLLCA 



ER 

FSV 

PI 

DY 

PW 



BE INITIALIZED AS YMAX(I) = M ^ X ( 1 , | Y ( 1 , I ) I) 

A VECTOR OF LENGTH MY 

A VECTOR OF LENGTH NY 

A VECTOR V LFNGTH N = NY 

A VECTOR jP LENGTH N = NY 

AN ARRAY IN WHICH THE J 

IN SUBROUTINE JACMAT WILL BE STORED. SIZE WHICH 
MUST 8£ ALLOWED IS DETERMINED BY THF STCRAGE TECH- 
NIQUE USED FOR IT, BUT NOFMALIY WON'T BE MORE THAN 
^ 2«M LOCATIONS, THE LA'^’^rp a^S'N BEING REQUIRED 
BY THE LINEAR EQUATION SOLVER. 



^ ML 
4- NL 

MA'^PIX CC-^PUTEC 



LCA 
LCA 
LDA 
LDA 
LCA 
LDA 
LCA 
LCA 
LCA 
LCA 
LCA 
LDA 



OIMENSTZN Y(7,l)» YUl), 5AVE(7,l), YMAX(l). ER(U, YLSV(l), FKDLOA 

1. PERT<A,3), CCJF(21), ESV(l), DY<1), PW<1), SAV(l). A(29i LDA 

EQUIVALENCE (A(8),BND), (A(9),BR), <A(10),r), ( A < 1 1) , ED W si ) , LCA 

l(A(12)f ENQl), ( A( l3) ,FNQ2 ) , < A ( 14 ) , ENP3 ) , (A<15),EPS), ( A ( 16 ) , EUP ) LD A 

2, < A< 17) ,HNtW ) , (Aa8),PEPSH), ( A ( 19 ) , I DOUB ) , ( A ( 20 ) , I WE VAL ) , LCA 

3<A<21)fK), (A(22) tLCOPYU , ( A ( 2 3 ) , LCOP V Y ) , C A ( 24 ) , MA X CL R ) , LCA 

4<A(25),M1), (A(26),NL), (A(27),NQ), (A<28),KSI, (A(29),Nvy) LCA 

LCA 

LDfl 

LDA 

TF5TING and lCA 
Dions. LDA 
LCA 

LCA 

»,l.,l.,.23, LCA 

LCA 

ICA 

LDA 

THE STIFFLY LCA 
MACHINE LCA 

LCA 
LDA 
LCA 
LCA 
LDA 
LCt 
LDA 
lD A 



TFE COEFFICIENTS IN THE PERT ARRAY ARE USED FjR tRPCR 
CHANGING STFPSIZE AND NEED TO BE ACCURATt. TC CNLY A F 

CAT A PERT/4.,9.,16.,25.,36.,49.,9.,16.,25.,:-c.,49.,64 
12. 7889E-2, 1.7C569E-3,6.83929E-5/ 

THE ENTRIES IN THE CCF ARRAY ARF THE COEFFICIENTS FCR 
STABLE METHCDS USED IN THIS PROGRAM AND ARE TO 5c THE 
PRECISION EQUIVALENTS OF THE FOLLOWING CGNSIANTS. 

-1 

-3/2 , -1/2 
- 11/6 , -1 , - 1/6 
-25/12 , -35/24 , -5/12 , -1/24 
-137/60 , -15/8 , -17/2‘t , -1/8 , -1/120 
-147/60, -2C3/90 , -49/48 , -35/144 , -7/240 , -1/720 



no 

eao 

ESO 
900 
910 
S2G 
930 
S40 
9 50 
960 
970 
9 80 
990 
ICOO 
1010 
1020 
1033 
1040 
1C5G 
1060 
1C70 
1080 
1090 
1100 
1110 
1120 
1130 
1140 
1150 
1 160 
1170 
1180 
1 190 
1200 
1210 
1220 
1230 
1240 
1250 
1260 
1270 
1280 
129C 
1300 
1210 
1320 
1330 
1340 
1350 
1360 
13 70 
1580 
1390 
1400 
lAlO 
1A20 
1 A3C 
1440 
1450 
1460 
1470 
1A80 
1490 
1500 
1510 
1520 
1530 
1540 
1550 
1560 
1570 
1580 
1590 
)6C0 
1610 
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C LDi 1620 

Q LDA 1620 

CiTA C0F/-1. ,-l . 5, -.5 t - 1. 8 33533, - 1.,-. 1666 667,-2. 0833 23 »-l. 458333, LC A 1640 

1- . 41666 6 7, -.04166 6 67, -2. 28 3 3 33, -1.8 75 ,- . 708233 3, - . 12 5 ,-.C0S 23 323 3 , LCA 16 50 

2- 2. 45, -2. 2 5 55 56, -1.0 206 33, -.243 0556, -.02 9 1666 7, -.00 13 83 8 39/ LD.4 16 60 

IF (JSTA'STI 100,110,150 LDA 1670 

LOA 1680 

IF THIS 15 A RESTART ENTRY, RESTORE Y A.'JO VL FROM THE S4VF ANC LCA 1690 

YLSV ARRAYS, ViHEP E THEY HERE SAVED BY A PREVIOUS CALL TC LOASAV. LCA 1700 

LOA 1710 

100 CALL COPYZ ( Y ,S A VE , LOOP YY ) LCA 1720 

CALL COPYZ ( YL, YLSV,LCOPYL ) LCA 172C 

GC TO 150 LDA 1740 

LOA 1750 

IF THIS IS THE FIRST CALL, INITIALIZE YW.AX, SCALE CERIVATlVcS, ANDLDA 1760 

INITIALIZE INDICATORS AND SET ORDER TO ONE. LOA 1770 

FOR OOUfiLE PRECIS ION, SET LCOPYL = 14+NY AND LC3PYL = 2*NL IF LOA 1780 

SLBRQUTINE CCPYZ IS IK SINGLE PRECISION. LCA 179C 

LDA 1800 

no KL - N-KY LOA 1610 

LCOPYY = 7*KY LCA 1820 

LCOPYL = NL LCA 1830 

R1 = PIKO(R,KY) LOA 1840 

EPS = SCRT ( FLOAT! HI n*RKSEPS LOA 1850 

MAXCEP = HIK0(MAXCR,6) LCA I860 

IF (IPPT.LE.C) GO TO 120 LOA 1670 

PRINT 3, N,KL,RMSEPS,TEND,H LCA 1880 

PRINT 4 LOA 1890 

120 KS = 0 LOA 1900 

NW = 0 LOA 1910 

LCA 1920 

CO 130 J-1,KY LOA 1930 

YHAX(J) - 4HAX1I1.,ABS(Y(1, J) )) LOA 1940 

120 Y(2,J) = Y(2,J)*H LDA 1950 

LDA 1960 

KC = 1 LOA 1970 

8R = 1. LOA 1980 

ASSIGN 190 TC IRET LCA 1990 

LOA 2000 

SET COEFFICIENTS FOR THE ORDER CURRENTLY BEING USED. LCA 2010 

E IS A TEST FOR ERRORS OF THE CURRENT ORDER NO LCA 2020 

tUP IS TC TEST FOR INCREASING THE ORDER, ECWN FOR DECREASING THE LCA 2030 

CRCER. LCA 2040 

LDA 2050 

140 K = NCT lKQ-1 )/2 LOA 2060 

CALL COPYZ (A(2),C0F(Ktl),N0) LDA 2070 

K = NC+1 LCA 2080 

IC0U8 = NC LCA 2090 

ENCl - .5/NC LCA 2100 

EKC2 * .5/K LCA 2110 

EKC3 = .5/(KC+2) LCA 2120 

PEPSH = EPS*»2 LOA 2130 

E = PERT (KC,1 )4PFPSH LCA 2140 

cUP = PPRT(KC,2»*PEPSH LCA 2150 

EOWN = PERT(KC,3)4PEPSH LDA 2160 

BNC =• ( PPS*ENQ3)T*2 LCA 2170 

IREVAL = 1 LOA 2180 

GC TO IRET, (190,200,490,570) LCA 2190 

150 IF (H.FC.HNEW) GO TO 190 LCA 2200 

LCA 2210 

IF CALLER HAS CHANGED H, RESCALE DERIVATIVES TO REFLECT THAT HNFW LCA 2220 

HAS USED CN THE LAST CALL. LCA 2230 

LDA 2240 

R = H/HKEH LDA 2250 

assign 190 TO IRET LCA 2260 

GC TO 61C LCA 2270 

LCA 2280 

SET JS"^ART TC NQ, THE CURRENT ORDER OF THE HETHTL), BSF3P= EXIT, LDA 2290 

AND SAVE THE CURRENT STEPSIZE IN HNFW. LCA 2300 

LDA 2310 

160 JSTART = NQ LOA 2220 

HNEW = H LCA 2330 

RETURN LCA 2340 

170 KS = NS+1 LCA 2250 

IF (IPRT.LE.O) GO TO 180 I CA 2260 



26 



ooooooo oooooooooo r»o o O OOOOO oo^ ooo ooo 



2CC 



PRINT DATA IF DESIRED BY USER 




LDA 


PRINT 1, NS,NW,NQ,H.T,(Y(ljI ) ,I=1,NY) 
IF (NL.GT.O) PRINT 2, (YL ( I ) , i^l,NL) 
CCNTINL-F 

IF (KFLAG. LT.O) GG TO 16H 
I F ( T.GE.TcNC ) GG TO 160 




LCA 

LDA 

LCA 

LCA 

LDA 


TAKE ANCTHtP STEP IF T < TEND 




LCA 


JSTART = 1 




ICA 


SAVE DATA FCP TRIAL WITH A SMALLER TIMESTEP IF THIS 


STEP FAILS 


LCA 


CALL COPYZ (SAVE, Y,LCGPYY) 
CALL COPYZ (YLSV, YL,LCOPYL ) 
RACUM = 1. 

KFLAG = 1 
HCLD = H 
NCOLD = NQ 
TCLD = T 
T = T+H 
HINV = l./H 




LDA 
LCA 
LCA 
LDA 
LCA 
LCA 
LDA 
LCA 
LCA 
— 1 r A 


CC.VIPUTP FSeClCftO VALUES PY SFFeCTTVPLY MULTIPLYING 
VECTCR PY PASCAL TPIANGL5 MA“PIX 


DERIVATIVE 


LCA 

LOA 


DC 210 J=2,K 
J3 = K+J-1 

DO 210 J1=J,K 
J2 = J3-J1 




LDA 

LCA 

LCA 

LCA 

LCA 

LDA 



OC 210 1=1, NV 

210 Y4J2,I) = Y(J2,n+Y( J2 + 1,I ) 



DC 220 1 = 1 , KY 
220 FR (I ) = C. 



DC 270 L=l,2 

C&lL DJPFUN ( Y, YL ,T,HINV,0Y) 
IF (I^PViSL.LT.l) GC TO 230 



IDA 
LCA 
LCA 
LOA 
LCA 
LDA 
LCA 
LCA 

LDA 

DC JO TC 7HPEE: CCfi'^ECTCk ITERATIONS. CGNVERGcNCE IS :BTAIN'FC WHENLDA 

CHANGES APE LESS THAN 6ND ^HICH IS DtP‘.= N0ENT ON THE EPkCF “^EST LCA 
CONSTANT, THE SJH OF COR'^ECTIONS 15 ACCUMULATED IN Efi(I). IT IS LCA 
EQUAL T: THE K-TH CERIYATIVE OF Y TIMES H*=*^K /( K-F AC TO R I AL *A ( K )) , LDi 
AND THUS IS PPOPOPT ION?L TO THr ACTUAL EFR0P5 TO THE L0r,-5T PGV^F^ LCA 
PF H PRESENT, WHICH IS H^-^'K. LCA 

lda 

! 

LCA 
LCA 
LDA 
-LCA 
LCA 
LCA 

LCA 

L C A 

LCA 
-LCA 
lCA 
LCA 
LCA 
LCA 
LCA 
LCA 
LCA 
LCA 
LOA 
LCA 
LCA 
LDA 
LCA 
LOA 



IF THERE HAS BEEN A CHANGE CF ORDER OR THERE HAS BEEN TRjUPLE 
WITH CONVERGENCE, PW IS R E -b V ALUATE’j PRIOR TO STARTING th- 
CCRPECTCR ITERATTCN. IWEVAL IS THEN SET TC -1 AS AN INCICATGR 
THAT IT HAS BEEN CONf=. NEWPW IS SET ^':NZER^. Tj iNClCATE TO 
^LBROUTINE KUITSL THAT A NEW ?W HAS BEEN PPCVlD'^O. 



CALL JACRAT ( Y , Yl , T , H i NV , A ( 2 ) , N , NY , E PS , DY , F 1 , PW ) 
KFLAG = 1 
IWEVAL = -1 
NW = NW+1 
N E W P W = 1 

230 CALL NLITSL ( PW , DY , FI , N,r4Y , E PS t YM AX ,NEwP w , Kpn cr ) 
IF (KRRET.NE.O) GC TO bOO 
IF (NL.Lc.O ) • GO TO 250 

CO 240 1^1, NL 

240 YL (I) = YL ( I )-Fl( I + NY ) 



25C CONTINUE 
CEL = 0. 



2 3 70 
2300 
23RC 
2400 
2410 
2420 
2420 
2^40 
2450 
2460 
2470 
2480 
24S0 
2500 
2510 
2520 
2530 
2540 
2550 
2560 
2570 
2580 
2590 
2 6 DC 
2tlC 
2620 
2630 
2640 
2650 
2660 
2670 
2680 
269C 
2700 
2710 
2720 
2730 

2 7 4C 
27CC 
2760 
2770 
2?B0 
2790 
2800 
2810 
2820 
2830 
2840 
2850 
2860 
2870 
2680 
2690 
2900 
2910 
2920 
2930 
2940 
2950 
2960 
2970 
2980 
2990 
3C0G 
3010 
3020 
3C30 
5040 
3050 
2060 
3070 
3080 
3090 

3 100 
3110 
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260 



270 



DO 260 1=1, N\ 

Y( 1,1) = Y 1, I)-F1( I ) 

Y(2,I ) = Y(2,I)>A(2)«F1( I ) 

ER(1) = ER(I)+F1( I) 

DEL = DEL+( FI (I )/AMAXl( YMAX( I ),ABS( Y( 1, I ) ) ) )^*2 
CONTINUE 

IF (L.GE.2) BR = AMAXl ( .9* BR , DE L/DEL 1 ) 

DELI = CE‘ 

IF (AMIM 
CONTINUE 



DELI = CEL 

IF (AMIM(0EL,BR^DEL*2.).LE.BND) GO TO 330 



280 

290 

200 

310 



T = TOLC 

IF (IWEVAL) 280,300,290 

IF (H.LE.HMIN^l. 00001) GO TO 310 

RACUM = RACING. 25 

CONTINUE 

GC TC 560 

KFLAG = -3 



RESTORE Y ANC YL AFTER CONVERGENCE FAILURE 



320 

330 

340 



call COPYZ (Y,S AVE,LC0PYY) 
CALL COPYZ (YL, YLSV,LC0PYL ) 
H = HOLC 
NC = NQCLD 
GO TO 170 



THE CORRECTCR CONVERGED, SO NC^^ THE ERROR TEST IS MADE. 



LDA 3120 
LDA 3130 
LD& 3140 
LDA 3150 
LDA 3160 
LDA 3170 
LDA 3180 
LDA 3190 
LDA 3200 
LDA 3210 
LDA 3220 
LDA 3220 
LDA 3240 
-LDA 3250 

TFE CCRRECTIQR ITERATION FAILED TO CONVERGE IN 3 TRIES. VARIOUS L CA 3260 
POSSIBILITI ES ARE CHECKED FOR. IF H IS ALREADY HMIN AND PW HAS LCA 3270 
ALREADY BEEN RE-EVALUATED. A NO CONVERGENCE EXIT IS TAKEN. LDA 3280 
CTHERWISE TFE MATRIX PW IS RE-EVALUATED AND/OR (IN THAT ORDER) THELDA 3290 
STEP IS REDUCED TO TRY AND GET CONVERGENCE. LDA 3300 

-LDA 3310 
LDA 2220 
LDA 3230 
LDA 3340 
LDA 3250 
LDA 3360 
LDA 3370 
LDA 3380 
-LDA 3590 
LDA 3400 
-LDA 3410 
LDA 3420 
LCA 3430 
LDA 3440 
LOA 3450 
LC.A 3460 
-LDA 3470 
LDA 3480 
-LDA 3490 
LDA 3500 
LDA 3510 
LCA 3520 
LCA 3530 
LOA 3540 
LDA 3550 
LOA 3560 
LOA 3570 
-LDA 3580 
LDA 3590 
LDA 3600 
LDA 3610 
LDA 3620 
LDA 3630 
LOA 3640 
-LDA 3650 
LDA 366C 
LDA 3670 
LDA 3680 
LDA 3690 
LDA 3700 
LDA 3710 
LDA 3720 
LDA 3730 
LDA 3740 
LCA 3750 
LDA 3760 
LDA 3770 
-LDA 3780 
LDA 3790 
LDA 3800 
LDA 3810 
LDA 3820 
-LDA 3830 
LDA 3840 
LDA 3850 
LDA 3860 



C = 0. 

DC 340 1=1, M 

YH = AMAXl ( AES(Y( 1, I ) ) ,YMAX( I ) ) 
D = D*( EP( I )/YM)**2 

IWEVAL = 0 

IF (C.GT.E) GO TO 380 



TFE ERRCR TEST IS OKAY, SO THE STEP IS ACCEPTED. IF lOOUR 
NCW BECOMES NEGATIVE, A TEST IS MADE TO SEE IF THE STEP SIZE 
CAN BE INCREASED AT THIS ORDER OR ONE HIGHER CR ONE LCWER. 

THE CHANGE IS MADE ONLY IE THE STEP CAN BE INCREASED BY AT 
LEAST 101. ICOUB IS SET TO NC TO PREVENT FURTHER TESTING 
FOR A while. if no CHANGE IS MADE, I DOUB IS SET TO 9. 

IF (K.LT.3) GO TO 360 

DC 350 J=3,K 

DO 350 1=1, NY 

350 Y(J,I) = Y(J,I)-^A(J)^ER(I) 

26C KFLAG = 1 

ICOUB = IDCLB-1 
IF (IDOLB) 410,370,510 
370 CALL COPYZ (ESV,ER,M1) 

GC TO 510 

THE FRRCR TEST FAILED. IF JSTART = 0, THE DcPIVATIVcS IN THE 
SAVE ARRAY ARE UPDATED. TESTS ARE THEN MADE ^0 FIX THE STEPSIZE 
AND PERHAPS REDUCE THE ORDER. AFTER RESTORING AND SCALING THE 
Y VARIABLES, THE STEP IS RETRIED. 



380 IF ( JSTART. GT.O) GO TO 400 
DO 390 1=1, NY 
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390 S^VE(2tI ) = \(2,I ) 



400 

410 

420 

4 30 
440 

45C 

460 

470 



KFLAG ^ KFLAG-2 
IF (H.LE.HMIN) GO TO 550 
T = TOLC 

IF (KFLAG.LE.-5 ) GO TO 530 
PR2 = ( C/E) ^♦ENQ2*1.2 
L * 0 

IF (NQ.LE. I ) GO TC 430 
C = 0. 

CC 420 J=lfM 

Y^^ = AMAXK AeS(Y( I, J ) ) , YMAX(J ) ) 

C = 0+( Y<K,J)/YM)**2 

PRl = (C/ECUK)*’>6NQ1*1.3 
IF (PR1.GE.PP2) GO TO 430 
PR2 * PRl 
L -1 

IF (KFLAG. LT.C. OR. NQ. GE. MAXOER ) GO TO 450 
C = 0 

DO 440 J=lfM 

YP = AMAXl (Ae$(Y( 1, J) ),Y^*AX( J)) 

C = 0+( (Eft( J)-ESV( J) )/YM)«’»'2 

PRl = ( C/EUF )*«ENC3*1.4 
IF (PR1.GE.FR2) GO TO 450 
PP2 = PRl 
L = 1 

R = l./ANAXl(PR2f 1.5-5) 

IF (KFLAG. LT.O. OR. R. GE. 1. 1 ) 



ICOUB = 9 
GO TO 510 
NEhC = NC+L 
K = NEWC-H 

IF (NEV^C.LE.NQ) GO TO 480 
R1 = A( KEWQ )/FLOAT(NEWQ) 

CC 470 J=ltNY 
Y(K,J) = ER(J)*R1 



GO TO 460 



480 CONTINUE 



IF THF STEP WAS OKAY, SCALE THE Y VARIABLES IN ACCORDANCE 
WITH THE NEW VALUE OF H. IF KFLAG < 0, HOWEVER, USE THE 
SAVED VALUES (IN SAVE AND YLSV). IN EITHER CASE, IF THE QRJER 
HAS CHANGED IT 1$ NECESSARY TO FIX CERTAIN PARAMETERS BY CALLING 
THE PROGRAM SEGMENT AT STATEMENT NUMBER 140. 

ICDUB = NC 

IF (N5wC.5Q.NC) GC TO 490 
NC - NEWC 

ASSIGN 490 TC IRET 
GO T<^ 140 

490 IF (KFLAG. GT.O) GO TO 500 
RACUM = PACUM^R 
GC TO 560 

50C R = AMAX1(AMIN1(HMAX/H,R) ,HMIN/H) 

H = H*R 
UEVAL = 1 
ASSIGN 510 TC IRET 
GC TO 610 

510 CC 520 1=1, M 

520 YMAX(I) = AHAXl (ABS( Y( 1,1 ) ) ,YMAX(I ) ) 

GC TO 170 

THE ERROR TEST HAS NOW FAILED THREE TIMES, SC THE DERIVATIVES ARE 
IN BAD SHAPE. RETURN TO FIRST ORDER METHOD AND TRY AGAIN. fF 
CCURSEf IF NC = 1 ALREADY, THEN THERE IS NC HOPE AND WE EXIT WITH 
KFLAG = -4. 

530 IF (NQ.EC.l) GO TO 540 



LCA 

LCA 

LDA 

LCA 

LDA 

LDA 

LCA 

LCA 

LDA 

LCA 

LDA 

LCA 

LCA 

LDA 

LCA 

LCA 

LDA 

LCA 

LDA 

LDA 

LCA 

LDA 

LCA 

LDA 

LCA 

LCA 

LDA 

LCA 

LCA 

LDA 

LCA 

LOA 

LDA 

LCA 

LDA 

LCA 

LOA 

LDA 

LCA 

LCA 

LCA 

LOA 

LDA 

-LDA 

LOA 

LDA 

LCA 

LCA 

LDA 

-LUA 

LCA 

LCA 

LCA 

LCA 

LDA 

LCA 

LDA 

LCA 

LCA 

LCA 

LCA 

LCA 

LDA 

LOA 

LCA 

LDA 

LDA 

LDA 

-LDA 

LCA 

LDA 

LDA 

LCA 

-LDA 

LDA 



3 8 70 
3880 
2890 
3900 
3910 
3920 
3920 
2940 
2950 
3960 
29 70 
3980 
5990 
4000 
4010 
4020 
4030 
4040 
4050 
4060 
4070 
4080 
4090 

4 100 
4110 
4120 
4130 
4 140 
4150 
4160 
4 170 
4180 
4190 
4200 
4 210 
4220 
9230 
4240 
4250 
4260 
4270 
4280 
4290 
4 300 
4 310 
4 3 20 
4230 
4340 
4350 
4360 
4270 
4380 
4390 
4400 
4410 
4420 
4430 
4440 
4 450 
9460 
4470 
4480 
4490 
4 500 
4 5 1C 
4520 
4530 
4540 
^550 
4560 
4570 
4580 
4 590 
4600 
4610 
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NC = 1 
ICOUB = 1 

ASSIGN 570 TC IRET 
GC TO lAO 
540 NCOLD = 1 
KFLAG = -4 
GO TO 220 
550 KFLAG = -1 
GC TC 170 



THIS SeCTirN RESTORES THE SAVED VALUES OF Y AND YL , SCALING ^ 
Y DERIVATIVES AS NECESSARY, AND THEN REiURKS TO THE PREDICTOR 



560 H = HOLC’frRACLM 

H = AMAXKHMN, ANIN1(H,HMAX) ) 
570 RACUM = h/hCLD 
R1 = 1. 

CO 580 J=2*K 
R1 = Rl’i'RACLM 



58C 



DC 580 I=1,KY 

Y( J,n = SAVE(J,I )=^R1 



CG 5S0 1=1, NY 
590 Y(l,n = SAVF(1,I ) 

CALL COPYZ (YL, YLSV,LCGPYU 
IREVAL = 1 
GC TO 20C 
6CC KFLAG = -5 
GC TO 16C 



THIS SECTION SCALES THE Y DERIVATIVES BY R**J 



6 10 R1 = 1. 

CC 620 J=2,K 
PI - R1«P 

CC 620 1=1, NY 
620 Y( J,n = Y ( J, n*Ri 

GC TO IRET, (190,510) 



THIS SECTION ALLOWS FOR RESTARTS AFTER -OLVING ANOTHER PROPLE 
HAVING TeP^«I^ATEO THE CURRENT COMPUTER kUN. SUBROUTINE LDA SA 
SAVES THE NECESSARY VALUES WHICH ARE INTERNAL TO LCASLB. FOR 
CCU6LE PRECISION. WITH COPYZ IN SINGLE PRECISION. THE NUMBER 
LCCATIONS TC PE SAVED J^ND RESTORED, LCOPYS AND LlOPYR. MUST fc 
SET TC 56. 

IT IS ASSUMED THAT IN ADDITION TO THE V4RIAPLES IN THE ARRAY 
SAVED BY CALLING LDASAV, THE USER ALSO SAVES THE ARRAYS SAVE, 
YLSV, YHAX, ESV, AND PW. 

TC RESTART THE USER FIRST CALLS LDARST TC RESTORE THE VALUES 
BY LDASAV, THEN RE-ENTERS LDASUB WITH JSTAPT < C, AND WITH TH 
CTHER PARAMETERS THE SAME AS RETURNED FROM THE LAST ENTRY Tj 
LCASUB, PARTICULARLY THOSE ARRAYS MENTIONED ABOVE. 

ENTRY LCASAY(SAV) 

LCOPYS = 29 

CALL COPYZ (SAV, A , LCOPYS » 

CALL COPYZ (SAVE, Y,LCOPYY ) 

CALL COPYZ (YLSV, YL,LCOPYU 
RETURN 

ENTRY LCARST(SAV) 

LCCPYR = 29 

CALL COPYZ (A,SAV,LCOPYR) 

RETURN 



LDA 
LCA 
LDA 
LCA 
LDA 
LDA 
LDA 
I CA 
LDA 
-LDA 
HE LDA 
LTDPLDA 

LCA 

LCA 

LDA 

LCA 

LDA 

LCA 

LCA 

LCA 

LDA 

LDA 

LDA 

LCA 

LDA 

LDA 

LCA 

LDA 

LLA 

LDA 

LDA 

LDA 

LDA 

LDA 

LCA 

LDA 

LCA 

LDA 

LDA 

LCA 

LDA 

LDA 

LDA 

LDA 

LCA 

LDA 

ORLDA 
LDA 
LDA 
LCA 
LDA 
LDA 
LDA 
LCA 
IDA 
LDA 
SAVECLDA 
LDA 
LDA 
LDA 
-LCA 
LDA 
LDA 
LDA 
LDA 
LCA 
LCA 
LDA 
LDA 
LDA 
LDA 
LDA 
LDA 
-LDA 
LDA 



M, 

V 

OF 

E 



462C 
4620 
4640 
4650 
4660 
467C 
468C 
4690 
4700 
4710 
4720 
4 7 30 

4 7 40 
4750 
4760 
4770 
4780 
4790 
4800 
4810 
4820 
4830 
4840 
4850 
4860 
4870 
4880 
4690 
4900 
4910 
4920 
4930 
494C 
4950 
4960 
4970 
4980 
4990 
5000 
5010 
5020 
5030 
5G40 
5050 
5060 
5070 
5080 
5090 
5100 
5110 
5120 
5130 
5140 
5150 

5 160 
5170 
5180 
5190 
5200 
5210 
5220 
5230 
5 240 
5250 
526C 
5270 
5280 
5290 
5200 
5310 
5320 
5230 
5340 
5350 
5260 
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1 FCf5f^Ar (2I5,I2» lo2E10.2,7El4.6/( 32X,7E14.6) i 

2 FORMAT C-2X ,1F7E14.6) 

3 FORMAT l*l N ='.I3,» NL =M3,‘ FM^EPS =*,lP69.2t» 

1 it9.2,' F =SE9.2//) 

4 FORMAT (• NS NW Q H'»8X,*T 

PNC 



Ti.N0 



•t8X,»Yil,«) AMD VL(^)’//) 



LCA 3270 
LCA 5280 
LCA 5290 
LCA 5400 
LQA 5^10 
LCA 5^20 
LCA 5430 



SL8R0UTIN0 CCPYZ(StY,L) 
DIMENSION S(l),Y(lJ 



This SUBROUTINE COPIES THE ARRAY Y, OP LENGTH L, INTO TPP AnFAY 



IF (L.Lc.C )RtTURN 
CG 100 J=ltL 
IOC S(J) = Y(J) 
RETURN 
PNC 



COP 
COP 
---CC-P 
rcp 
S COP 
CCF 
■--OOP 

COP 

COP 

COP 

CCF 



10 

20 

30 

40 

50 

60 

70 

80 

90 

LOO 

110 

120 



SLBROUTINE CERVAL ( Y t YL t T t N t N Y , U t KF R£T ) 



DER 
DGR 
CEO 



C 

c 



THIS SUBRCUTINF CALCULATES THE INTIAL VALUES uH THE CFRTVAT*VPS 
IN THE GENERAL CASE. IT IS WRITTEN SO THAT IT SHOULD WORK IF THP CEP 
FIRST NY EOLATIONS ALL INVOLVE DERIVATIVES. 1“ ATTEMPTS TO SOLVE CEP 
THE FIRST NY EQUATIONS USING NEWTON'S METHOD, BUT SINCE IT TRIES CEB 
TC EVALUATE CF/DY' BY CALLING JACMAT IN SUCH A WAY AS TO MAKL the OER 
CF/DY TERM INSIGNIFICANT, IT IS POSSIBLE THAT IT MAY FAIL FOP tha-^OEP 
REASON. IT MAY FAIL FOR OTHER REASONS, AS WELL. IF IT DOES FAIL CEP 
THE USER CAN SUPPLY HIS uWN VERSION OF DERVAL, OP MODIFY TH.5 CEP 

ROUTINE IN SUITABLE FASHION. THIS ROUTINE ASSUMES THAT VALUES CF CEP 
THE linear VARIABLES HAVE BEEN SUPPLIEC PREVIOUSLY. IF THOSE CFR 

MUST BP SCLVEC FOR SIMULTANEOUSLY WITH THE CEPiVATIVtS, THF USER DER 
VLST SUPPLY HIS OWN VERSION OF DERVAL. DFP 

OPR 
DEF 
DEP 
CER 
OtR 
CER 
DER 
DEO 
Dtp 
CER 



THE CALLING SEQUENCE FOR THIS SUBROUTINE IS 

CALL D'^PVAL(Y,YL,T,N,N'r,W,KCRET J 

WHERE THE PARAMETERS ARE DEFINED AS FOLLOWS 



YL 

T 

N 

NY 

W 

KERET 



SAME AS IN LDASUB AND SOESOL. Y(l,n CONTAINS THE 
INITIAL VALUES OF THE DEPENDENT VARIABLES. THF 
VALUES OF THE DERIVATIVES ARE OETURNED IN Y ( 2 , . 

SAME AS IN LDASUB ANC SOESOL. THE INITIAL VALJCS OF CEP 
THE LINEAR VARIABLES MUST BE SUPPLIED TO THIS yrRSIClNuEF 
initial TIME DEP 

SAME AS IN LDASUB, TOTAL NUMBER OF VARMELES CER 

SAME AS IN LDASUB, NUMBER OF CIFFFRENTIAL EQUATIONS DEP 
AND NONLINEAR VARIABLES DEP 

SCRATCH ARRAY W FROM THE CALLING SEQUENCE OF SDESCL. CtP 
THIS CAN BE USED AS NEEDED IN this SUBROUTINE. DCP 

RETURN INDCATOP DER 

=0 NORMAL PETURN DER 

=1 ERROR PETURN DER 

DER 

CER 

DEB 
CcR 
DER 
DPR 
Dtp 
DER 
DEP 
PER 
DER 
DtR 
DEP 
DER 
DEP 
DcR 
DER 



DIMENSION Y(7,l), YL(1), W(l) 

CC 100 1=1, NY 

W(2<‘N + I) = AMAX1(ABS<Y(1,I) 1,1.) 
ICO Y(3,I ) = 0. 

HINV = 16. **20 
KERET = C 
EPS2 ^ NY/1.E8 
EPS = SCRT(EPS2) 

DC 140 IT = iaO 

DC 110 1=1, NY 
llO Y(2,I ) = Y(2,I)/HINV 



10 
20 
30 
40 
50 
60 
70 
80 
90 
100 
lie 
120 
12G 
140 
150 
160 
170 
180 
190 
200 
210 
220 
230 
240 
250 
260 
27C 
280 
290 
20C 
310 
220 
3 30 
240 
250 
2 60 
270 
28C 
390 
^00 
410 
420 
420 
^40 
450 
^60 
470 
^80 
^90 
500 
510 
520 
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120 



130 



140 



150 

160 



170 



Ct^LL DIFFUN ( Y , YL , T , HI NV , W ) 

Ci^LL J£C^^AT <Y,YL.T,HINV,-1.,NY,NY,EPS,W,W(N^1) ,W(3’<‘N + 1) ) 
NEUPVi - 1 

CC 120 I^lfNY 
U( I) = V^d )’*HINV 

CALL NUIT5L (W(3 + N-H ) » W , W ( N + 1 ) , NYt NY , C PS , U< 2<^N+1) » NF W Fl^ , K RE 7 ) 
IF (KPPT.NF.O) GO TG 170 
ER = 0. 

DC 130 1=^1, NY 

Y(3,n Y(3,I)-W(N*I) 

ER = EP+(W(N+I)/AMAXl(ABS(Y(3fI ))»1.) )*<=2 

IF (EP.LT.EPS2) GO TO 150 
CONTINUE 

GC TO 17C 

CC 160 1=1, NY 
Y(2,I) = YO,I) 

RETURN 
KERET = 1 
RETURN 
PNC 

SLBF3UT INE JACMAT ( Y , YL , T , HI NV , A2 ,N , N Y , E P S , DY , F 1 , P ) 



SLBPCUTINE JACMAT IS (USUALLY) SUPPLIED BY THE USER. ITS P 
IS TO EVALUATE THE J MATRIX NEEDED WHEN THE CORRECTCR EQU 
IS solved by NEWTCN'S METHOD. THIS VERSION APPROXIMATES 
J BY NU^^ERICAL CIFFERENCING AND USES FULL STORAGE MOLE 
IN AN NXN Matrix. 



URPCS5 

mTIGN 



JACMAI 



:alculates the matrix 



J = 



DF A2 DF 

CY H DY» 

THE CALLING SEQUENCE FOR THIS SUBROUTINE IS 
CALL JACMAT ( Y.YL,TiHlNV,A2 ,EPS,N,NY,DY,F1,PW) 

where the parameters are defined as FOLLOWS. 

Y - SAME AS IN LDASUB AND IN SDESOL. ON INPUT TO FFIS 

SUBROUTINE THE ARRAY CONTAINS CURRENT VALUES or THE 
DEPENDENT VARIABLES AND THEIR (SCALED) DERIVATIVES. 

YL - SAME AS IN LDASUB AND IN SDESCl. ON INPUT TO THIS 

SUBROUTINE THE ARRAY CONTAINS CURRENT VALUES OF THE 
LINEAR VARIABLES. 

T - CURRENT TIME 

HINV - 1/H , WHERE H IS THE CURRENT STEPSIZE 

A2 - A(2) FROM LDASUB. 

N - SAME AS IN LDASUB, TOTAL NUMBER OF VARIABLES 

NY - SAME AS IN LDASUB, NUMBER OF DIFFERENTIAL EOUhTIONS 

AND NCNLINFAR VARIABLES 

EPS - L2 ERROR CONSTANT USED IN LDASUB; 

OY - ARRAY OF FUNCTION VALUES A* CURRENT VALUES OF THE 

VARIABLES, INPU"»* TO JACMAT- 

FI - SCRATCH ARRAY OF N LOCATIONS WHICH CAN BE USED EY 

THIS SUBROUTINE IN ANY WAY NEtCED. 

PW - J MATRIX, OR APPROXIMATION, CALCULATED IN JACMAT AN 

RETURNED TO CALLING PROGRAM. THIS MATRIX IS USED IN 
SUBROUTINE NUITSL AND STORAGE MODE MUST AGREE EETWEE 
THE TWO SUBROUTINES. 

DIMENSION DY(1), Y(7,l), YL(1), FlU), °W(1) 



DER 


5 30 


DEP 


540 


DER 


550 


OER 


560 


DER 


570 


HER 


580 


DER 


5<50 


DEP 


600 


HER 


610 


DSR 


620 


DEP 


630 


DER 


640 


DcP 


650 


DER 


660 


DEP 


670 


DEP 


680 


DEP 


690 


DEP 


700 


DER 


710 


DER 


720 


DER 


730 


DEP 


740 


DER 


7 50 


DER 


760 


DER 


770 


DER 


780 


DEP 


790 


DER 


800 


JAC 


10 


-JAC 


2? 


JAC 


3C 


JAC 


40 


JAC 


50 


JAC 


60 


JAC 


70 


JAC 


80 


JAC 


90 


JAC 


100 


JAC 


no 


JAC 


120 


JAC 


130 


JAC 


140 


JAC 


150 


JAC 


160 


JAC 


170 


JAC 


180 


JAC 


190 


JAC 


200 


JAC 


210 


JAC 


220 


JAC 


230 


JAC 


240 


JAC 


250 


JAC 


260 


JAC 


270 


JAC 


280 


JAC 


290 


JAC 


200 


JAC 


210 


JAC 


220 


JAC 


330 


JAC 


340 


JAC 


250 


JAC 


260 


JAC 


2 70 


JAC 


330 


DJAC 


290 


JAC 


400 


NJAC 


410 


JAC 


420 


JAC 


430 


-JAC 


440 


JAC 


450 
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M = N-NV 
= N^N 
C 

CC 100 1=1, NN 
100 = 0 . 
c 
c 

CC 12C J=1,NY 
F = Y( 1,J J 
E = Y(2,J) 

= EPS+^^^^X1(EPS,ABS(F),ABS(E) ) 
Y(l, J) = Y( 1,J)+R 
Y(2, J) = Y(2, J)-A2^R 
call OIFFUN (Y, YL ,T,HINV,F1) 

C 

CC 110 1=1, N 

110 PVy(! + ( J-1)*N) = ( Fin )-0Y( I n/R 
C 

Y(2, J) = F 
12C Y(1,J) = F 
C 

TF (M. EC.O ) GO VO 150 

c 

CC 140 J = 1,M 
F ^ YL(J) 

R = FPS^AMAX1(EPS,A0S(F) ) 

YL(J) = YL(JJ+R 

CALL DIFFUN ( Y, YL , T , HI N V , F 1 » 

C 

CC 130 1 = 1, N 

130 Pk(I+(J4NY-l)«N) = (Fl(I)-JY(I J )/R 
C 

140 YL(J) = F 
C 

150 CCNTIN'UE 
RETURN 
EMC 



J AC 
JAC 
JAC 
lAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 



460 
470 
480 
490 
500 
510 
520 
530 
540 
550 
560 
570 
580 
590 
600 
6 10 
62C 
630 
64C 
650 
660 
670 
680 
690 
700 
710 
720 
730 
740 
750 
760 
770 
78C 
790 
6CC 
810 
820 



SUBROUTINE NUITSL ( PW, OY , F 1 , N ,NY , EPS , Yr'AX , NEWPW , KR ET ) 

THE FUPFCSE OF THIS SUBROUTINE IS TO SOLVE A 
LINEAR SYSTEM Of EQUATIONS FOP THE NEWTON ITERATES WHEN THE 



NUI 

-NUI 

NUI 

NUI 



CORRECTuP EQUATION IS BEING SOLVED. UPON ENTRY TO THIS SUB aC UT I N ENU I 



THE SYSTEM CF EQUATION'S TO BE SOLVED IS 
J IS STORED IN PW UPON ENTRY 
W I S RETURNED IN FI 
-F IS STCRED IN DY UPON ENTRY 



J W = 



WHFRE 



NUI 

NUI 

NU! 

NUI 

NUI 



THIS SUBROUTINE IS GENERALLY SUPPLIED BY THE USER, ALTHOUGH THERE NUI 
ARE SOME STANCARO FORMS AVAILABLE. FOR EXAMPLE, iHlS VERSION 
ASSUMES THAT PW IS STORED IN FULL STOPAC-E MODE IN AN NXN MATRIX. 

IF NEWPW = 1, AN LU DEC GM PO 5 ! T I CN IS DONE, NEWPW IS SET *^0 ZERO 
A^,C FDRWARC AND BACKWARD SUBSTITUTION FOR THE SOLUTION IS DONE. 

IF NEWPW == 0, ONLY FORWARD AND BACKWARD SUBSTITuTTCN FOR THE 
SOLUTICN IS NECESSARY. 



NUI 
NLI 
NU! 
NUI 
NUI 
NUI 
NUI 

NCTE THAT THIS VERSION OF NUITSL REQUIRES THAT PW HAVE N^ + 2 + 2-r« NUI 
LCCATICNS SINCE 2 *N LOCATIONS ARE US^D BY THE IMSL LINEAR EQLATIONNUI 
SCLVEP. NLI 

NUI 

NCTE that THE PARAMETERS EPS AND YMAX ARE USEFUL IF aN ITERmTIvE NUI 
METHOD IS USED TO SOLVE THE SYSTEM OF EQUATIONS. NUI 

NUI 

THE CALLING SEQUENCE FOR THIS SL'BROU*riNE IS NUI 

NUI 

CALL NUITSL (PW,DY ,F 1 , N , NY , E P S , YMA X, NE WP W , KP ET ) NUl 

NUI 

WHERE the parameters ARE DEFINED AS FOLLOWS. NLI 

NUI 

PW - THE J MATRIX CALCULATED IN SUBRCUTINE JACMAT NUI 

DY - THE RIGHT HAND SIDE OF THE LINEAR SYSTEM TO BE SuLVEDNUl 

FI - THE SOLUTION IS RETURNED IN THE ARRAY FI NUI 

N - SAME AS IN LDASUB, TOTAL NUMBER OF VARlABLtS NUI 

NY - SAME AS IN LDASUB, NUMBER CF C I FFERE N*^ I A L FQJATICNS NUI 
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100 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 
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270 
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320 

^30 

340 

350 
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EPS 

YMAX 

NEWPH 



KRET 



AND NONLINEAR VARIABLES 
L2 ERROR CONSTANT USED IN wDASUB 

MAXIMUM VALUES OF Y(1,I) SEEN UP TO THE CURRENT TIME 
INDICATES WHETHER A NEW J MATRIX HAS BEEN COMPUTED 
=1 INDICATES A NEW J MATRIX HAS BEEN COMPUTED 
SINCE THE LAST ENTRY TO NUITSL. NEWPW 
SHOULD BE SET TO ZERO IF SOME PREPROCESSING 
SUCH AS LU DECOMPOSITION MUST BE DONE TN A 
NEW J MATRIX. 

»0 INDICATES THE J MATRIX IS THE SAME AS WHEN 
NUITSL WAS LAST ENTERED 
RETURN INDICATOR 

=0 NORMAL RETURN 

=^1 ERROR RETURN. SOLUTION OF EQUATIONS CCULD 
NOT BE OBTAINED. 



CIMENSICN PW<1), DY(l)f Fid), YMAX(l) 

NL = N-NY 

IF (NEWPW. EC. 0) GO TO 100 
NEWPW = 0 
NN * N* + 2«H 
NNN * NN+N 

CALL LUCATF (PW • P W , N.N . 0, D 1 , D2» PW (NN ) » PW ( NNN ) , F I , I ER ) 
IF ( lER.EO.O) GO TO io6 
KPET » 1 
RETURN 

IOC CALL LUELMF ( PW , D Y , PW ( NN) , N , N , F 1) 

KRET = 0 

RETURN 

EKC 



NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

'NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 

NUI 



SUBROUTINE Cl FFUN ( Y , YLtT , HINV f DY ) 

ilBROUTINE cTffuN IS SUPPLlio BY THE USER. ITS PURPOSE IS TC 
EVALUATE THE FUNCTIONS AT CURRENT VALUES OF THE VARIABLES. 

THE CALLING SEQUENCE FOR THIS SUBROUTINE IS 

CALL DIFFUN(V,YLiTfHINVf DY) 

WHERE THE PARAMETERS ARE DEFINED AS FOLLOWS. 

Y 

YL 
T 

HINV 
DY 



SAME AS IN LDASUB AND SDESQL. ON INPUT TO THIS 
SUBROUTINE THE ARRAY CONTAINS CURRENT VALUES OF THE 
DEPENDENT VARIABLES AND THEIR (SCALED) DERIVATIVES. 
SAME AS IN LDASUB AND SDESOL. ON INPUT TO THIS 
SUBROUTINE THE ARRAY CONTAINS CURRENT VALUES OF THE 
LINEAR VARIABLES. 

CURRENT TIME 

1/H . WHERE H IS THE CURRENT STEPSIZE 
RETURNED ARRAY OF FUNCTION VALUES. 



DIMENSICN Y(7,1),YL(1),DY(1) 

DEFINE YOLR FUNCTION HERE 

RETURN 

EKC 



270 
380 
290 
AOO 
410 
420 
430 
440 
450 
460 
470 
480 
490 
500 
510 
520 
530 
540 
550 
560 
570 
580 
590 
600 
610 
620 
630 
640 
650 
• 660 
670 
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Appendix 2: Examples 



Example 1: This example is the problem proposed by Gear [3]. The 

system of equations is 



^ S 



+ (R-y^)^ + I 

j=l J J 



= 1, 2, 3, 4 



1 r 

where R = ^ L y • 
i=l ^ 



1 ^ 9 

s = i I (R-y.)^ , 

i=i 



and 



75 + 7i yj + 7i 76 = 0 



^^6* H ' ^1* ''l' ^ ■ ° 



'>1 - '>2 * H fe ‘ ° 



+ V2 + Sy^ 72 = 0 . 



In the above b-- = b^^ = b^« = b,, = 447.501 
11 22 33 44 



'>12 ■ -‘’34 - ‘21 = -‘43 - - “ 2.499 



‘l3 ■ -‘24 - ‘31 ' -‘42 - - “-“5 



‘l4 - ‘23 ‘41 -‘32 - • 



The initial conditions are 

Vj = 1 , 1-1, 2. 3, 4 . 

75-70-1 

Vi = - 2 , V 2 = - 3 . 

9F 

Note that a different version of DERVAL is necessary since — is singular. 

9y 
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CIMCNST :N Y(7,6I,YL(2),W(1^)0),3! (161 
/C6T/G(16) 

DATA GI/<»^7.501 -^52 . A9 9 , 7 . 499 , -52 . 5 01 1 -452 . 499 , 447 .5 Ji ,52, 

1 4 7.49 9,-4 7.499, 5 2. 501,4*+?. 501,452.49 9,-52.501,47.499,452.4 9^ 

2 447.501/ 

OAT A N, NY,M t ^,REPS ,HMAX, HMIN,H, T, TEND/8, 6 ,2,6, 1.6-2, 5.62,1.: 
1 l.E-4,0,, 1.F3/ 

CALL EPPScT (207,256,-1, 1 ) 

CALL EPPSET (208,256,-1,1 ) 

CG 8 1^1,16 

8 cm ^ GI (I ) 

DC 10 1=1,4 
1C Y( 1,1) = -1. 

Yd, 5) = 1. 

Y(l|6) = 1. 

YL( 1) = -2. 

YL(2) = -2. 

JSKF = 0 

call SOESGL (Y,YL,T,T6NO,NY,NL,'^, JSKF,6, l,H,HWliN,HviAX, ^EP'£, W) 
PRINT 6,J$K6 

6 FCRMAT( *0 JSKF =* ,14) 

STOP 

ENC 



SLBPGUTINE CIFFUN ( Y , YL , T , h I .nIV , DY ) 

CC?MjN /CAT/G(4,4) 

DATA TGLC/-12.459/ 

OINE.NSI > Y(7,l ), YL( 1),DY(1 ) 

IF (T.EO.TCLOiGC TG 10 
TNTERM = EXP(-T) 

TCLD = T 
10 CCNTINUE 

R = (Y(l,l) 4 Yd, 2) + Yd, 3) + Yd,4))/2. 

S = 0. 

DC 20 1=1,4 

20 $ = S (R - Y( 1, I) )*=J‘2/2. 

DC 30 1=1,4 

DY(I) = FlNV*Y(2,I) - S +(R - Y(1,I))«^2 
CC 25 J=l,4 

25 DY(I) = CY(I) ^ G(I,J)>^Y(1, J) 

3C CONTINUE 

CY(5) = *-INV*(Y(2,5) + Y ( 1 , 1 ) *Y ( 2 , 6 ) + Y ( 2 , 1 ) ^ Y( 1 , 6 ) ) 

CY(6) = 2.=«'Y(1,6) + Y(l,6)^>!«3 - Yd,l) * YL (1 ) - l.-T^TER^ 
CV(7) = YL(1) - YL(2) Y ( 1 , 1 ) ^ Y ( 1 , 6 ) 

CV(8) = YL(1) -I- YL(2) ^ 5.^Y(1,1) * Y(l,2) 

RETURN 

ENC 



SLBPGUTINE CERV AL ( Y. YL , T , N, N Y, ^ ,KEPET ) 
CIMENSICN Y(7,l ), YL( 1) ,!^(1 ) 

KERcT = 0 
DG 50 1=1, NY 
5C Y(2,I) = C. 

HINV = 1. 

CALL CIFFUN(Y,YL,T,HINV,W) 

CC 100 1=1, NY 
ICC Y(2, I) = -^d ) 

RETURN 

END 



5C1 , 
, 
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Example 2: This example is a small one, contrived to illustrate the 
possibility of derivatives entering in a nonlinear fashion. The equations 
are 



Yl - 98y^ + 98y2 = 0 
^2 ■ ® ^ Yi + 199y2 = 0 




The initial conditions are 



= Y2 = 1 = 2 . 

Note that we have supplied the explicit expression for the Jacobian. 
Either JACMAT or a modified version of DERVAL must be supplied as the 
numerical difference approximation to the Jacobian causes DERVAL to 
fail. 
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OIMfNSICN Y(7,2).YL(l),W(50) 

DATA N,NY,NLt^tRtPS,HMAX,HMINtH,T,TEN'C/3, 2, 1 , 2f 1. E- 5 , 25 . , 1 . S- 10, 
1 l»E“4fC»f5C»/ 

Y(l,l) r 1. 

Y( 1,2) = 1. 

YUl) = 2. 

JSKF = 0 

CALL SDESCL ( Y , Y L , T ,T END ,N Y , NL , M, JSKF, 6, 1 , h, HMIN, HM AX, REPS, W ) 
PRINT 6, JSKF 
5 FCRMAT( »0 JSKF ,14) 

STCP 

ENC 



SLBROUTINE CIFFUN ( Y , YL , f , HINV, 0 Y ) 

DIMENSION Y(7,1),YL( l),DY(l) 

DATA T0LC/-79.03/ 

IF(T.EQ.TOLC)GC TO 10 
TNTERM = EXP(“T ) 

TCLD = T 
10 CONTINUE 

DY(1) = Y( 2,1)’*'HINV - 98.*Y(1*1) + 99.*Y<1,2) 

DY(2) = (Y( 2,1)*HINV )<‘«2 ♦ Y(2,2)*HINV - (196. - ^ MT 6 RM )^Y ( 1 , 1 ) + 
1 199.*Y(1,2) 

DY(3) ^ VL( 1) - Y (1,1) - Y( 1,2) 

RETURN 

ENC 



SL6RDUTINE JACMAT ( Y , YL , T, HI NV, A2 tN, NY , E '^5 , CY, FI , P W ) 
CIMENSTCN Y(7,1),YL(1),F1(1 ) , D Y ( 1 ) , P ( N , 1 ) 

AH = A2«HINV 
DG 100 1 = 1, N 
CC 100 J=liN 
100 PVs( I , J) ^ 6. 

PU(1,1) 

Pl^ (1,2) 

PW(2,1) 

PW(2,2) 



-AH 

9S. 



- 98. 



IF(NL.LE.O)RETURN 



-2.*AH^Y( 2,1 )*HINV - 198. + EXP(-T) 
-ah * 199. 



P^i(2,l) 

PW(3,2) 

PM3,3) 

RETURN 

ENC 



1. 
1 . 

= 1 . 
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Example 3: This is another contrived example, this one to illustrate 

the use of one type of sparse matrix storage, along with the use of 
iteration to solve the equations (2.4). The system of equations is 



A y + B y = 0 , 



where 

I 7 0 0 -3 

2 8 0 0 

__ 0 0 10 
A = 

-3 0 0 5 

0-100 
^ 0 0 0 -2 



.3 0 0 .1 

13 0 0 

0 0 10 



B = 



10 0 0 20 

0 5 0 0 

^0 0 0 57 

The initial conditions are 



= 1000 
Y2 = 0 

Y3 = -25 
= 10 

75 = 0 

= -1000 



0 

0 

0 

0 

4 

0 

0 

0 

0 

0 

6 

0 



0 

0 

0 

0 

6 

-.2 

0 

0 

0 

0 

100 
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The matrix storage scheme used for A, B, and the Jacobian, 
since it has nonzero elements in the same positions as A and B , is 
that outlined by Gustavson [4]. Briefly, one stores a pointer array 
(here called JS) which indicates the initial position of new elements 
in two other arrays, one of which (here called JN) gives the column 
number of the element stored in the corresponding position of the 
coefficient arrays (here called A and B) . 

Thus, for the above problem the arrays stored are 



JS: 


1 4 


6 ,^^7^ 


9 


ii__ 


_13 


























JN: 


1 4 


6^2 






^ 4 




T~ 


■^6 4 


A : 


7 -3 


-1 8 


2 


1 


5 


-3 


4 -1 


6 -2 


B : 


.3 .1 


-.2 3 


1 


1 


20 


10 


6 5 


100 57 


The elements of the 


i— 


row 


are 


Stored beginning 


at location 


JS(i) 


of the 


array A , 


and in particular 


A(JS(i) + k) is the element 


in the 


1— row and JN(JS(1) 


+ k) 




column of A , 


for k = 0 , 1, . . 



JS(i+l) - JS(i) - 1 • For our purposes it is necessary to access the 
diagonal element easily, so we have required that the diagonal element 
be the first element stored for a given row. This means that 
JN(JS(i)) = i , i=l, ..., n . Note that JS(i) is the number of nonzero 
elements in rows 1 through i - 1 , and that JS(n+l) must be defined as 
the total number of nonzero elements. 

Problems similar to the above arise when the finite element method 
is used to discretize the space domain for time dependent partial differ- 
ential equations. Simple modifications to the subroutines given below 
should permit solution of large problems arising in that fashion. We 
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note, however, that it is not convenient to store synmetric matrices 
in this form unless all nonzero elements are stored. Storage of only 
the elements of the lower triangular matrix requires one to reference 
columns of the matrix, which are not readily accessible. Even if the 
entire matrix is stored, total storage requirements for matrices arising 
in finite element applications is still considerably less with this 
scheme than that required by symmetric band storage mode [5]. 
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ICO 

110 

12C 



3 

4 

5 

6 
7 



DINtNSI^N Y ( 7,6 I , w( 1261 *Y I ( 6 ) , 12 ) , 12 ) , JSC(7 ) » JNC( 12 ) 

C:MM3i\‘ /CAT4/-fl(12) ,B(12),N, JS(7) ,JN(12) 

J6,JN 

CATA T, TEND,H, JSKF / 0 . ,2 30 . , 1 . F- 5 , 0 / 

CATA JSC/1,A,6,7,c,11,13/ 

CATA JNC/ I, 3, 5, 2, 1,3,^, 1,5, 2, 6, 4/ 

DATA AP/7# ,6*,~2*/ 

DATA BC/.3,.l,-.2,3.,l.,l.,20.,10.,o.,5.,lC0.,57./ 

CATA YT/1000.,0.,-25. , 10^ , 0. ,“1000./ 

\ = 6 

N Fl = f\ 4 1 
JE = JSC(NFl) “ 1 
CC 100 1^1, N 
Y(1,I ) = Yt (I ) 

cc 110 1 = 1 , je 

A(I) = AC(I ) 

en) = ec(i) 

JN(I ) = JNC (I ) 

OZ 120 1=1, NPl 
JS (I ) = JSC (! ) 

PRINT 4,N, ( JS(I ) ,I=l,NPl) 

PRINT 5,(Ji\(I),I = l, JE) 
point 7,(AU),I=l,JE) 

PRINT 6j(e( I ) ,I =1, JE) 

CALL SDcSOL (Y,YL, f ,T6Nij,N ,0 ,N, JSKP, 6, l,H,l.F-6,l*i5.,l.E“4,W) 

PRINT 3,JSKF 

STCP 

cCRf''AT( ‘CRETLRN FRGM SCcSGL TH JSKF =»,I4) 

FORMAT!’ FOR THIS CASE \=*,I3//’ fHF JS ARP AY V/ ( I 2 T 10 ) ) 
FGRMAT(/’OTFF JN AP R A Y • / / ( I 2 1 10 ) ) 

FCRMAT(/*0THE 3 ARRAY’//! 12F10.2) ) 

FORMAT ! /’OThE A A RR AY • / / ! 1 2 F 10. 2 ) ) 

PNC 



SCBRjUTINc ClfcjN ! Y, YL,T, HINV,DY ) 

CCMMON /CATA/A!12) ,8!12),N, JS!7), JN!12) 

INTEGER»2 JS,JN 

CIMENSICN Y !7, 1 ) , YL! 1) ,0Y! I ) 

CC 400 1 = 1, N 
CY!I ) = C. 

JE = JS ! I ) 

JS = JS !I + 1 ) - 1 
DC 300 J=je,JE 

CY!I) = CY!l) 4 y !2, JN! J) )=t'A!J)’4'Hli>lV 4 3 ! J ) « V ! 1 , J N ! J ) ) 
200 CCNTINUE 
400 CCNTINLE 
RETURN 
PNC 



SL3RGUT INS NUITSL ! PW , DY , F 1 , N , N Y , EP S , Y Me, X , N W , K R ) 
CCMM3N /C AT A/ A! 12 ) ,8! 12 ) ,ND, JS !7) , JN! 12 ) 

INTEGER42 JS,JN 

CIMENSICN Pa! 1) ,CY! 1) , FI! 1 ), YMAX! 1) 

CATA CMEG,C^EG^1 /I. 05, .05/ 

KRET = 0 
SPSS = FFS^^2 
6F$A2 = EPPS’S*. 0001 
NCIT = N 
CC 100 I=1,NY 

ICO FI !I ) = CyI I ) /o^i JS! I ) ) 

CC 300 IT=1,NCIT 
RCh = 0. 

CF = 0. 

CC 200 1 = 1, NY 
je = JS !i ) 4 1 

JE = JS!I41) ^ 1 
FN = DY!I) 

IF! JB.GT.JE )GC Tj 18C 
CC 150 J=JE,JE 

150 FN = FN - PVs!J) =«‘F1! JN( J) ) 

130 FN = FN/FW!je-l) 

FN = “ F1!I)*GM8GM1 

ACH = FI !I ) - 
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CH = CH 4 ( ACri/Y^AX( I) 

RCH = PCH 4 ( ACH/AMAXUABS (fN),6P^ ) 
200 FUT) = FN 

IF (RCH.LT.EPSS) RETURN 
IF(CH.IT.£PSA2) RETURN 
300 CCNTINUE 
KRET = 1 
RETURN 
ENC 



3ieROUTIN6 JACMAT(Y. YUfT^HINV, A2»N,NY»EPS, CV,F1»PW ) 
CCMMGM /CATA/A(12) »6( 12 ) , NDUI^, J S ( 7 ) , JN ( 12 ) 

INTEGER*2 JStJN 

Clf'ENSICN Y(7,l)tYL( DfFK 1),CY< 1),PW(U 
AH - -A2*HINV 
JE = JS(N4U - 1 
CC 100 J=1»JE 

ICO PU(J) s ^H*A(J) + B(J) 

RETURN 

ENL 
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Example 4: This example arises from a nonlinear reactor dynamics 

problem where the finite element method is used to discretize the space 
domain. The resulting system of equations has the form 

• 

A y - B y + w(C y) y = 0 , 



where C is a matrix with three subscripts. The i — equation can be 
expressed as 



N 



N N 



I 

j=i 




In this example 
per row in A 




N = 28 
and B . 






+ w 



I 

3=1 



I 

k=l 



'ijk 



^3 



yv - 0 



, and there are at most seven nonzero elements 
The nonlinear term y . y, appears only if 

1 K 



both a^^ and are nonzero. Therefore a different type of sparse 

matrix storage is used for this problem. 

An array, K , dimensioned (28, 7) is used to store (for each row), 
the columns subscripts for the nonzero elements. For convenience in 
accessing the diagonal element, we require that K(i,l) = i . We can 
note this matrix is simply the connectivity matrix for the finite 
element grid. Then the nonzero elements of A and B , are stored in the 
corresponding portions of the arrays A and B respectively. If there 
are in fact less than seven nonzero coefficients in a row, the remaining 
K(i,j) are set to zero. 

The storage for C is somewhat more complicated. C is symmetric 
(invariant under any permutation of subscripts). The nonlinear term of 



th 

the i — equation was rewritten as 
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N N _ N N _ 

w I ^ \ik ^ ^ ^iik ^k ’ 

j=l k=l j=l k=j ^ 



where 



ijk 



‘'ijk 

^Ijk ^ikj 



k < j 
k = j 
. k > j 



The coefficients then stored in a (28,28) array C in the 

order the second and third subscripts are given here. 

(K(i,l), K(i,D), ..., (K(i,l), K(i,7)), (K(i,2), K(i,2)), (K(i,2), 

K(i,7)) (K(i,7), K(i,7)) . 

The equations can then be written in the form 



7 7 



- hj yKCi.:)’ * J, ^KCij) yu±, 



k) 



= 0 



where 



„ . k + UzlHiizJi . 

Jk 2 



Because of the large amount of data for this problem the input 
arrays K, A, B, and C are simply listed along with the programs for 
this example. 
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Clf-ENSICN Y(7,28) ,WS(560) 

CCMMON /CAT4/M28,7) tB(28,7),C(2e,20) ,N,NNZ,K(28,7) 

IMEG£«^2 K 

CATA TPNO,H^IN,HM/5X,EPS,ZOM£GA/ .1, l.E“12t .1 ,.01 t 650.90 3£-U/ 
C/ITA NN,NY,NL/28, 28tO/ 

^^Z = 7 

C/iLL EPRSET (207,256,-1,1 ) 

C4LL EPRSET(2C8,256,-1,1) 

N = NN 
N = N 

NE^C == NNZ+ (NNZ ♦ l)/2 
PRINT 10 

cc no m,NN 

PE4D 1, (K( 1 ,J) ,J = 1,NNZ) 
no PRINT llt(K(I,J) ,J = 1,NNZI 
PRINT 12 
CC 120 1=1, NN 
RE/iD 2, (A(t ,J) ,J=1,NNZ) 

V(1,I) = 0. 

120 PRINT 15,(A(I, J» , JxI.nNZ) 

Y(l,l) = 1.E16 

PRINT 13 

CC 130 !=1,NN 

RE/iC 2, (8(I,J),J = 1,NNZ) 

130 PRINT 15,<e(I, J) , J=1,NNZ) 

PRINT 14 
CC 140 1 = 1, NN 
RE^D 2, (C( I,J) , J=1,NENDI 
CC 135 J=1,NEN0 
135 C(I,J) = d(I,JI ^ZOMEGA 
140 PRINT 15,(C(I,J), J=1,NEND) 

JSKF = 0 
T = 0. 

H = HMIN^IOOO. 

CALL SDESOL(Y,YL,T,TEND,NY,NL,K, J5KF,6,l,H,h>«IN,H!^AX, FPS, ,>£ ) 

PRINT 3, JSKF 

STCP 

1 FCPMAT(16I5) 

2 FORMAT! 7E11. 4) 

3 FORMAT ( «0 JSKF = • , I3» 

1C FCRMAT( *1K ^RPAY« /) 

11 FORMAT (8X, 1418) 

12 FCRMAT(///»0A(I , J )•/) 

13 FORMAT!///* 06(1 ,J)*/I 

14 FCPMAT!///*CC!I,J)*/) 

15 FCRMAT!8X,1F7E16.6) 

END 



SLBROUTINE CIFFUN ! Y, YL,T,HINV,DY I 

CCMMON /CATA/A! 28,7) ,B!28, 7) ,d( 28,28) ,N,NNZ,K(28, 7) 

INTEGERY2 K 

CIMFNSKN Y!7,1),YL!1),DY(1) 

DO 400 1=1, N 
CY!I ) =0, 

CO 300 J=1,NNZ 

IF (K! I, J).LE.O)GC TO 310 

QY!I) = CY!I) ♦ Y!2,K(I,J) )«A(I,J)+HINV - B ( I , J ) * Y ! 1 , K ( I , J ) ) ♦ 
1 C!I,J)^Y(l,K!l,J))=)'Y!l,K!I,l)i 
30C CCNTINUt 
310 L = NNZ 

CC 360 J1=2,NNZ 
IF!K(I,J1).LE.0)GC TO 400 
DC 350 J2=J1,NNZ 
L = L + 1 

IF(K! I, J2).LE.0)G0 TO 350 

DY(I) = CY(I) + C!I,L)^Y!1,K(I,J1))>J‘Y!1,K!T,J2)) 

350 CCNTINUB 
360 CCNTINLH 
400 CONTINUE 
RETURN 
END 

SLBPOUTINE JACMAT !Y,YL, T, HI NV,A.2,N, NY , EPS, CY,F 1,PW I 
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COMMON /CAT4/A( 28t 7) tR(28, 7) 28,28) ,NO,NNZ,K (28, 7) 
INTEGFR^2 K,P 

DIMENSION Y(7,Uf YL(1),F1( I ) , D Y (1 ) , P W ( N Y , 1 ) 
CIMENSI1N P (7,7) 

DATA P/AS*0/, lNITP/0/ 

IF (INITP.EC.NNZ)GG TQ 99 
IMTP = NNZ 
CO 98 L=1,NNZ 
DC 98 M = LfNNZ 

98 P(L,M) = M + (L - 1)*(2*NNZ - L)/2 

99 CCNTINUE 

AH =-A2«hlNV 
CO 300 1=1, NY 
DC 300 J=1,NNZ 

PWd , J) = Ah*A( !, J) - &( I , J ) 

CC 100 1=1, J 

IF(K(I,L).L£.0)G0 TO 295 
ICO Pi^(I,J) = PVn(I,J) + C(T,P(L, J))*Y(1,K(I,D) 

CO 200 M=J,NNZ 
IF(K( I , M).LE.O)Gu TO 295 
200 PU(IfJ) = PW(I,J) C(I ,P( J,M) )«Y(1,K(I,M) ) 

295 CCNTINLE 
300 CCNTINUE 
RETURN 
END 



SUBROUTINE NUI TS L ( , DY , 1 , N , NY , E PS , Y MA X , N E WP W , KR E T ) 
DIMENSION PUlNYd ) ,0Y(1) ,F1(1 ),YMAX(1 ) 

CCMMCN /CATfl/A(28,7) ,8(28, 7) ,C (28,28) ,ND,NNZ,K(23, 7) 
CATA SPC,SPCM1/1 .05,.05/ 

INTEGEP«2 K 
KRET = 0 
EFSS = EFS'«'«2 
EPSA2 = EPSS'^.OOOl 
NCIT ^ N 

280 CC 281 1 = 1, NY 

281 Fl( I ) = CY( I )/PW( 1,1) 

CC 287 IT=1,N01T 

RCF = 0, 

CF = 0. 

CO 285 1=1, NY 
FN = DY(I) 

CO 28A J=2,NNZ 

IF(K( I, J).LE.O.OP.KC, J).GT.NY)GO TJ 28A 
FN = FN - PVsd, J)+F1(K(I, J) ) 

28A CONTINUE 

FN = FN/FU(I,1) 

FN = FN*SPb - SPDM1*F1( I ) 

ACH = FI (I ) - FN 
CH = CH 4 ( ACH/YMAX( I ) )*’!‘2 
PCH = RCF 4 (ACH/AMAX1(AB$(FN),EP$) )‘^‘’^2 
285 Fid) = FN 

IF(RCH.LT.EPSS )Cn TO 288 
IF(CH.LE.EPSA2)GG TO 288 

287 CCNTINUE 
KPET = 3 

288 CCNTINUE 
RETURN 
END 
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INPUT DATA FOR EXAMPLE 4 



16 2 
2 1 6 



4 


3 


e 9 


10 5 




















5 


4 


1C 21 






















6 


11 


7 2 


1 




















7 


2 


6 11 


12 13 




a 
















8 


3 


2 7 


13 9 




4 
















9 


4 


8 13 


14 15 




10 
















10 


5 


4 9 


15 22 




21 
















11 


6 


16 12 


7 




















12 


7 


11 16 


17 18 




13 
















13 


8 


7 12 


18 14 




9 
















14 


9 


13 16 


19 20 




15 
















15 


10 


9 14 


20 23 




22 
















16 


11 


28 17 


12 




















17 


12 


16 28 


27 18 




















18 


13 


12 17 


27 26 




19 14 
















19 


14 


18 26 


25 20 




















20 


15 


14 19 


25 24 




23 
















21 


5 


1C 22 






















22 


21 


10 15 


23 




















23 


22 


15 20 ■ 


24 




















24 


23 


20 25 






















25 


20 


19 26 


24 




















26 


19 


16 27 






















27 


18 


17 28 


26 




















28 


16 


27 17 






















8.3776E 


02 


8.3776E 02 


4.1888E 


02 


0.0 




0.0 




0.0 




0.0 




5.0265E 


03 


4.1888E 02 


2.09445 


03 


2.5133E 


03 


2.09445 


03 


4.18885 


02 


Q.O 




1.6755E 


03 


4.1888E 02 


1.6755E 


03 


4.188 8E 


02 


0.0 




0.0 




0.0 




5.0265E 


Q3 


4.1888E 02 


2.0944E 


03 


2.5133E 


03 


2.09^4E 


03 


4.1888E 


02 


0.0 




1.4661E 


03 


4.1888E 02 


1 . 466 IE 


03 


3.1416E 


02 


0.0 




0.0 




0.0 




1.C891E 


04 


2.9322E 03 


4.1888E 


03 


2.0944E 


03 


8.37765 


02 


0.0 




0.0 




2.8484E 


04 


2.5133E 03 


4.1888E 


03 


6.2832E 


03 


6.7C21E 


03 


6.2832E 


03 


4.1883E 


03 


2. 1782E 


04 


1.6755E 03 


2.0944E 


03 


4.1888E 


03 


5. 8643E 


03 


4.1688E 


C3 


2.0944E 


03 


2.64E4E 


04 


2.5133E 03 


4.1888E 


03 


6.2832E 


03 


6.7021E 


03 


6.2832E 


03 


4.18385 


03 


1.9059E 


OA 


1.46615 03 


2.0944E 


03 


4.1888E 


03 


5. 1313= 


03 


3.1416E 


03 


1.5708E 


03 


2.3457E 


04 


2.9322E 03 


5.0265E 


03 


8.3776E 


03 


6.2832E 


03 


o«o 




0.0 




5.3616E 


04 


6.7021E 03 


8.37765 


03 


1.0472E 


04 


1.C891E 


04 


1.0472E 


O-^ 


8.37765 


02 


4.6S14E 


04 


5.8643E 03 


6.2332E 


03 


8.3776E 


03 


1.0053E 


04 


8.3776E 


03 


6.2832E 


03 


5.3616E 


04 


6.7021E 03 


8.3776E 


03 


1.0472E 


04 


1.C891E 


04 


1.0^725 


04 


8.3776E 


03 


4.4663E 


04 


5.1313E 03 


6.2832E 


03 


8.37765 


03 


8.9535E 


03 


8.4823E 


0-» 


6.2046E 


03 


3.2515E 


04 


5.0265E 03 


5.18365 


03 


1.0812E 


04 


1. 0472E 


04 


0.0 




0.0 




5.62C8E 


04 


1.C891E 04 


1. C812E 


04 


1.1958E 


04 


1. 1958E 


04 


1.0812E 


04 


0.0 




8.0582E 


04 


1.0053E 04 


1.0472E 


04 


1.0812E 


04 


1.33125 


04 


1.3312F 


04 


2.12845 


04 


5.62C8E 


04 


1.08S1E 04 


1.C812E 


04 


1.1958E 


04 


1.1958E 


04 


1.0812E 


04 


0.0 




6.C750E 


04 


8.9535E 03 


1.0472E 


04 


1.0812E 


04 


9.2481E 


03 


1.0348E 


04 


9.8764E 


0*3 


3.7699E 


03 


3.14165 02 


1. 5708E 


03 


1.8850E 


03 














2.9531E 




1.6850E 03 


3.1416E 


03 


6.2046E 


03 


8.71796 


03 










6.1889E 


04 


8.7179E 03 


8.4823E 


03 


9. 6764E 


03 


1. 2468E 


04 










6.3303E 




1.2466E 04 


1.03f8E 


04 


8.8357E 


03 














5.7962E 


04 


9.2481E 03 


1.1958E 


04 


1.4726E 


04 


8986575 


03 










5.4719E 


5^ 


1.1958E 04 


1.3312E 


04 


1.7671E 


04 














9.4719E 


04 


1.3212E 04 


1.1958E 


04 


1.4726E 


04 


1.7671E 


04 










5.3014E 




5.1836E 03 


1.4726E 


04 


1.1958E 


04 














-1.5816E 


09 


1.17205 09 


1.0449E 


09 


0.0 




0.0 




0.0 




0.0 




-3.9622E 


2^ 


1. 04495 09 


6.3535E 


08 


4.43385 


09 


3.3535E 


08 


1.0^49E 


09 


0.0 




-3.1631E 


09 


1.0449E 09 


2.3440E 


09 


1.0449E 


09 


0.0 




0.0 




0.0 




-3.9E22E 


09 


1.0449E 09 


6.3535E 


08 


4.4338E 


09 


6. 3535E 


08 


1.0449E 


09 


0.0 




-4.3361E 


09 


1.0449E 09 


1.8355E 


09 


1.4879E 


09 


0.0 




0.0 




0.0 




-6.7925E 


09 


4.56C9E 09 


6.7778E 


09 


6.35355 


08 


1. 1720E 


09 


0.0 




0.0 




-1.5222E 


10 


4.4338E 09 


6.7778E 


09 


1.9061E 


09 


1.1212E 


10 


1.9061E 


09 


6.7776E 


09 


-1.3585E 


1C 


2.3440E 09 


6.3535E 


08 


6.7778E 


09 


9.12185 


09 


6.77785 


09 


6.3535E 


08 


-1.5223E 


10 


4.4338E 09 


6.7778E 


09 


1.9061E 


09 


1.1212E 


10 


1.90615 


09 


6.7778E 


09 


-2.4104E 


10 


1.83555 09 


6.3535E 


08 


6.7778E 


09 


7.3355E 


09 


8.4446E 


09- 


-6.0318E 


08 


-1.3995E 


10 


4.5609E 09 


7.9498E 


09 


1.3556E 


10 


1.9061E 


09 


0.0 




0.0 




-2.9627E 


12 


1.1212E 10 


1.35565 


10 


3.1768E 


09 


1.7989E 


10 


3.17665 


09 


1.5556E 


iO 


-2.7989E 


10 


9.1218E 09 


1.9061E 


09 


1.3556E 


10 


1.5900E 


10 


1.3556E 


10 


I.9061E 


09 


-2.9627E 


10 


1.1212E 10 


1.3556E 


10 


3.1768E 


09 


1.7989E 


10 


3.17685 


09 


1.3556E 


10 


-4.8828E 


10 


7.3355E 09 


1.9061E 


09 


1.3556E 


1C 


1.0212E 


10 


1.1621E 


10 


2.0406E 


09 


-3.5210E 


10 


7.9498E 09 


1.3692E 


10 


1.6043E 


10 


3.1766E 


09 


0.0 




0.0 




-8.2408E 


10 


1.7989E 10 


1.6043E 


10 


-3.6945E 


08 


2.4060E 


10 


1.3103E 


1C 


0.0 





-7.2321E 


10 


l.?9C06 


10 


3. 17686 


09 


1.3103E 


10 


1. 1476E 


10 


1.14766 


10 


1.6280E 


10 


-8.2^C8E 


10 


1.7989E 


10 


1.3103E 


10 


2. 4060c 


10- 


•3.6945E 


08 


1.6C43L 


10 


0.0 




-e.711SE 


10 


1.0212E 


10 


3. 1768E 


09 


1.60436 


10 


2.4798E 


10 


3.46583 


09 


1.33966 


1C 


-8.2637E 


09 


1 .48796 


09- 


6.03186 


08 


2. 89536 


09 














-A.7C95E 


10 


2, 89523 


09 


8 . 4446E 


09 


2.04036 


0 9 


6.4C05E 


08 










-1.0136E 


11 


6.40C56 


08 


1.1621F 


10 


1.3398E 


10 


4.68226 


09 










-1.3280c 


11 


4.68226 


09 


3.4658E 


09 
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